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GEMACS  ENGINEERING  MANUAL 


A.  INTRODUCTION 

The  function  of  this  manual  is  to  discuss  the  physical  and  mathema- 
tical methods  used  in  the  GEMACS  (General  Electromagnetic  Model  for  the 
Analysis  of  Complex  Systems)  code  to  obtain  the  results  desired.  There 
are  five  basic  steps  involved  in  obtaining  the  observables  used  to  evalu- 
ate a particular  configuration.  These  are: 

(1)  Geometric  Modeling 

(2)  Structure  Excitation 

(3)  Interaction  Computation 
(I4)  Numerical  Solution 

(5)  Observable  Computation 

Section  B discusses  each  of  these  areas  as  they  pertain  to  the  GEMACS  code. 
Section  C is  a discussion  of  the  BMI  (Banded  Matrix  Iteration)  solution 
technique  developed  as  a result  of  this  project.  Section  D is  a summary 
of  results  obtained  with  BMI  for  various  geometries. 

Much  of  the  information  in  section  B is  presented  in  the  AMP 
manuals.  Much  of  the  AMP  discussion  is  directly  applicable  since  the 
basis  formalism  and  some  of  the  code  itself  is  identical  to  the  AMP 
code.  The  primary  differences  are  in  the  geometry  and  excitation 
processes  and  in  the  use  of  the  BMI  technique  to  obtain  a solution. 

B.  GEMACS  FUNCTIONS 

When  the  GEMACS  program  was  first  structured,  very  little  was  known  of 
the  solution  techniques  to  be  implemented.  Therefore,  the  basic  functions 
needed  to  perform  a generic  electromagnetic  analysis  were  identified. 

The  input  requirements  for  each  function  were  identified  and  this  specified 
the  output  requirements  of  the  logically  preceding  function.  In  order 
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to  assure  flexibility  and  modularity,  the  outputs  of  each  function  are 
symbolically  identified  data  sets.  These  data  sets,  along  with  other 
parameters  are  used  as  inputs  to  subsequent  functions. 

Thus,  each  function  will  operate  on  previously  defined  data  or  user 
supplied  inputs  to  generate  another  data  set  for  subsequent  use. 

1 . Geometric  Modeling 

Geometric  modeling  is  used  to  convey  structure  geometry  to  the 
GEMACS  code.  Since  GEMACS  uses  the  EFIE  (Electric  Field  Integral  Equation), 
the  actual  structure  is  represented  by  a series  of  wire  segments  connected 
in  a way  that  approximates  the  actual  surface.  In  the  usual  solution 
technique,  the  segments  are  given  a number  and  the  number  assigned  to 
one  segment  is  totally  independent  of  other  segment  numbers.  This  is 
also  true  in  GEMACS  when  using  the  full  interaction  matrix  in  the  solution; 
however,  when  using  the  BMI  technique,  it  becomes  important  to  number 
the  segments  in  such  a way  that  segments  which  are  electrically  close 
have  numbers  which  are  also  close.  This  is  because,  for  any  given  row 
or  column  of  the  interaction  matrix,  a segment's  position  is  determined 
by  its  number.  Thus,  all  interactions  with  segment  1 will  be  in  row  1 
and  column  1.  Much  study  on  numbering  has  been  completed  and  is  summa- 
rized in  section  C.  In  order  to  accommodate  the  need  for  specific  numbering 
sequence,  the  RN  (Renumber)  command  has  been  included  in  the  geometry 
processor.  With  this  command,  the  user  may  enter  the  model  in  the  most 
convenient  manner  and  then  specify  the  desired  numbering  sequence.  When 
reducing  the  EFIE  to  a set  of  simultaneous  linear  equations,  two  assump- 
tions are  made: 

(1)  The  wire  segments  are  of  sufficiently  small  radii  that  circum- 
ferential currents  may  be  ignored. 

(2)  The  current  on  any  segment  may  be  approximated  by  some  function 
called  the  basis  or  expansion  function. 

These  two  assumptions  have  several  implications  depending  on  the  structure 
being  modeled  and  the  expansion  function  used.  Assumption  (1)  restricts 
the  wire  radius,  R,  to  something  less  than  a wavelength  (X)  . Typically, 
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R < A/10  is  quoted  as  acceptable.  If  the  expansion  function  is  to 


represent  the  actual  curren  over  the  segment,  then  the  segment  length, 
L,  must,  in  general,  also  be  less  than  A.  The  basis  function  used  in 
GEMACS  is  the  sine  + cosine  + pulse  expansion  with  collocation  and  it 
has  been  observed  that  L = A/4  is  sufficient  where  the  current  does  not 
vary  rapidly  and  A/10  is  adequate  where  rapid  variations  occur.  A 
general  rule  of  thumb  for  modeling  structures  is  that  the  total  area  of 
the  wire  segments  should  approximate  the  surface  area  of  the  structure. 
In  addition,  for  the  expansion  function  used,  a ratio  L/R  ^ 5 has  been 
found  to  give  consistent  agreement  with  currents  obtained  from  analytic 
solutions.  Using  these  rules,  the  number  of  segments  (N)  for  a given 
frequency  may  be  approximately  determined  from: 


N 

Y.  2nR. L.  = A 


where  A is  the  surface  area  to  be  modeled.  With  N = 800,  an  area  of 
2 

approximately  10  X can  be  represented  quite  accurately.  If  less  confi- 
dence in  the  observables  is  acceptable,  A may  be  much  larger  than 
indicated.  The  reader  is  cautioned  that  the  applicability  of  these 
rules  to  codes  employing  other  methods  is  largely  unknown  at  this  time. 

A consistent  method  of  modeling  has  not  been  established  except  for  the 
present  code. 

Once  the  structure  has  been  defined  to  the  GEMACS  code,  several 
operations  take  place  before  the  next  command  is  executed.  The  geometry 
data  are  scanned  to  find  all  segments  connected  to  each  other.  During 
this  operation,  segments  which  have  identical  end  points  are  found 
and  all  but  one  of  these  segments  are  given  a segment  and  tag  number  of 
zero.  In  this  way,  segments  in  planes  of  symmetry  and  on  axes  of  rota- 
tion are  allowed.  The  criterion  for  segments  being  connected  is  that  the 
end  point  separation  is  less  than  the  roundoff  error  of  the  host  com- 
puter. This  number  is  internally  computed  based  on  the  number  of  bits 
used  for  the  mantissa  of  a floating  point  number.  During  this  time,  any 
segment  which  terminates  on  the  XY  plane  is  identified  as  being  connected 
to  a ground  plane  if  one  is  subsequently  specified.  If  no  connections 
are  found,  this  is  also  noted  for  user  convenience.  Once  all  of 
the  junctions  are  found,  a list  is  constructed  which  identifies  the  next 
segment  connected  to  either  end  of  each  segment.  For  example,  if  end  1 
of  segments  1,  3>  5>  20,  and  75  are  connected,  then  data  identifying  the 
next  segments  would  be  3,  5.  20,  75,  and  1.  This  forms  a circular 
linked  list  identifying  all  segments  connected  to  either  end  of  a given 
segment.  This  information  is  listed  in  the  geometry  output  for  both  end 
1 and  end  2 of  each  segment.  A negative  number  implies  end  1 of  the 
identified  segment  while  a positive  number  identifies  end  2.  Once  all 
junctions  are  identified  and  the  junction  linked  list  constructed,  the 
geometry  data  internal  format  is  changed  from  end  point  data  to  centerpoint 
coordinates,  segment  length,  and  direction  cosine  format.  When  used, 
dimensioned  geometry  data  are  scaled  to  wavelength  for  computational 
ease. 
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Structure  Excitation 


There  are  two  methods  of  structure  excitation  available  in  the 
GEMACS  code.  These  are  voltage  excitation  for  antennas  and  a wave  excita- 
tion for  scatterers.  These  types  may  be  superimposed  for  composite  excita- 
tions . 

a . Antenna  Source  Model 

The  normal  derivation  of  the  currents  on  a body  involves 
writing  the  electric  field  as 

E(X)  = E mC  (X)  + E^X)  (Eq.  2) 

where  E (X)  is  the  total  electric  field  at  X,  E 1 nc  (X)  is  the  incident 
electric  field,  and  E (X)  is  the  electric  field  due  to  some  volume 
distribution  of  currents.  By  imposing  the  boundary  condition  at  a 
finite  number  of  points  on  the  surface  of  a perfect  conductor  and  speci- 
fying the  incident  field,  we  solve  for  the  currents  induced  on  the  body. 

(Eq.  3) 


(Eq.  4) 


n x E 


( Tine  ^ TS  \ 

= n x I E +E  1=0 


therefore, 


~ “inc  * ~S 
n x E = -n  x E 


In  the  thin-wire  approximation,  this  condition  is  enforced 
on  the  component  of  the  field  in  the  direction  of  the  wire  since  the 
assumption  that  azimuthal  currents  are  negligible  forces  the  0 components 
to  zero,  i .e. , 


l 


— i nc 


(Eq-  5) 
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where  11  is  a unit  vector  in  the  direction  of  the  wire.  This  condition 
is  enforced  at  the  centerpoint  of  each  segment  and  no  information  about 
E is  known  between  these  points.  This  is  equivalent  to  throwing  all  of 
the  wires  away  since  wire  parameters  have  no  effect  on  the  boundary 
conditions,  but  only  on  the  current  distribution  used  to  satisfy  the 
boundary  condition. 

Now,  the  only  way  to  specify  a boundary  point  to  the  code 

is  via  a wire  segment  and  thus  an  antenna  source  must  be  specified  as  a 

A — 

wire  segment.  However,  this  segment  is  totally  artificial  in  that  nXE 
along  this  segment  is  not  zero  since  V=J*  E • dl  is  the  voltage  driving 
the  antenna.  Since  the  boundary  condition  is  only  satisfied  at  the 
midpoint  and  the  structure  of  the  excitation  field  is  of  second  order,  a 
uniform  excitation  gives  reasonably  good  results  as  long  as  the  excited 
segment  length  is  much  less  than  X.  Since  the  total  field  at  the  midpoint 
has  been  specified,  equation  2 can  be  used  by  again  invoking  the  thin- 
wire  approximation  and  taking  only  the  tangential  field. 


‘ ' Eant  * ‘ (Ei"C  + £S) 

(Ei"C  - Eant)-  ^ • E$ 


■«.  • E' 


(Eq.  6) 


with  E = 0 at  those  points  not  driven  as  antenna  sources, 

ant 

The  solution  of  equation  6 is  a set  of  currents  at  the 
midpoints  of  each  segment  and  the  antenna  input  admittance  can  be  approxi- 
mated by: 


Y*  V 


(Eq.  7) 


The  shortcomings  of  this  type  of  model  are  presented  in 
the  AMP  Engineering  Manual  and  its  references  (reference  lb). 
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b . Field  Excitation 

The  incident  electric  field  (E)  on  a structure  may  be 
either  a polarized  or  nonpolarized  field  from  a plane  or  spherical  wave. 
Multiple  wave  excitations  may  be  combined  in  any  order  desired. 

In  general,  the  field  at  the  observation  point  r is  given 
by 


¥(r)  = ¥.  e"jk'  " +E  e"jkr 

! K 


(Eq.  8) 


where  E.  is  the  source  or  direct  illumination  field  and  E is  the  reflec- 
I K 

ted  field.  If  there  is  no  ground  plane,  E = 0.  There  are  two  differ- 
ences  between  a spherical  and  plane  wave.  The  spherical  wave  has  a 1/R 
amplitude  dependence  and  the  wave  vector  k is  always  oriented  in  the 
direction  of  the  field  point  from  the  source  point.  For  a plane  wave, 
the  amplitude  and  k are  constant. 

In  general,  the  incident  or  source  field  is: 


E , = E + j EP 


(Eq.  9) 


where  EP  is  the  pc'?*i:atio'  component  if  present.  _EP  is  determined 

— — a I EP  I 

from  the  vector  relationship  E x EP  = k and  ECC  = 4=--1'  where  ECC  is  the 
user  supplied  eccentricity.  |E] 

The  reflected  field  is  given  by 


F-R  -R11  F11  + Ru  F11  " r-l  Ej- 


(Eq.  10) 


Rj j is  the  in-plane  reflection  coefficient 
RA  is  the  out-of-plane  reflection  coefficient 


Ejj  is  the  in-plane  tangential  component 
Ej'j'  is  the  in-plane  normal  component 
Ea  is  the  out-of-plane  component 


* 


k 


Let  p be  a unit  vector  normal  to  the  plane  of  incidence.  Then,  with  k 

i 

from  source  to  specular  point  and  kr  from  specular  point  to  the  field 
point,  we  have 

p = k.  x kr  (unit  vector  1 to  plane  of  incidence) 

E = E,,  + Ex  (Eq.  11) 


Ex  = (E  • p)  p 

|-k  (E  + j EP  ) + k (E  + j EP  )1 
= L V x J x x v J y J 


(kx2  + kV2> 


["ky  ' * kx  J ] 
(Eq.  12) 


E11 

= E - Ej. 

(Eq. 

13) 

n 

Ei, 

= (En  • Z)  Z 

(Eq. 

14) 

En 

= E - E P 
11  bll 

(Eq. 

15) 

The  reflection  coefficients  R^ ^ and  are  the  modified 
Fresnel  coefficients.  This  method  is  discussed  at  some  length  in  the 
AMP  Engineering  Manual  (reference  16).  Portions  of  the  manual  are  repro- 
duced as  appendix  A.  The  reflection  coefficients  are: 


(Eq.  16) 
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L 


\ 


I 

t 


R 


1 


e 


E 


cos  0 - 


cos  0 + 


(Eq.  17) 


(Eq.  13) 


where  0 is  the  angle  of  incidence  measured  from  z and  Ej  and  are  the 

dielectric  constant  and  conductivity  of  the  ground  plane.  The  free  space 

permittivity  is  c . 

o 

3 • Interaction  Computation 

The  interaction  or  structure  matrix  used  in  the  GEMACS  code  is 
formally  identical  to  that  found  in  the  AMP  code.  Minor  changes  in  the 
data  structure  have  been  made;  however,  given  identical  geometries,  one 
obtains  identical  elements  in  the  structure  matrix.  The  derivation 
involved  in  obtaining  the  elements  of  the  structure  matrix  is  presented 
in  chapter  II  of  the  AMP  Engineering  Manual  portions  of  which  are  repro- 
duced here  with  equation  numbers  changed  to  be  consistent  with  numbering 
in  this  manua 1 . 

a.  Integral  Equation  Formulation 

The  electric  field  E due  to  a volume  current  distribution 
J is  written  by  means  of  the  Green's  dyadic  as 


E (ro)  - 


j u)yQ  J (r)  • 


G (r,  rQ)  d V (Eq.  19) 


where  rQ  and  r are  the  observation  and  source  points,  respectively,  and 
the  Green's  dyadic  is  expressed  in  the  usual  notation  as 


G (r,  rQ)  = -(1/^tt)  •[ I + (l/k2)w]  g 

where 

g = exp  (-jk  |7  - 7 | ) / |7  - 7 j 

o o 
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1 


and  I is  the  unit  second-rank  tensor.  The  suppressed  time  variation 
is  exp  (jwt)  with  cj  the  radian  frequency.  The  plane  wave  propagation 
constant  is  k,  and  is  related  to  eQ  and  the  permittivity  and  per- 
meability of  free  space,  respectively,  and  u>  by 


l 

I Where  the  current  distribution  is  limited  to  the  surfaces 

of  a perfectly  conducting  body,  equation  19  becomes 

E (?o)  = JJ  jwuo  Js  (7)  • G (7,  7Q)  dA  (Eq.  20) 

I S 

with  J the  surface  current  density.  If  this  surface  current  is  induced 
7 . —I 

by  an  incident  electric  field  E , then  an  integral  equation  for  the 
^ unknown  surface  current  can  be  obtained  from  equation  20  and  the 

boundary  condition  that 

n(rQ)  x [lfS  (7q)  + I1  (7q)J  = 0 (Eq.  21) 

' - _ - -S 

where  n (r  ) is  the  unit  normal  vector  at  r and  E is  the  scattered 

° ° - S 

field  due  to  the  secondary  current  distribution.  Equating  E of 

equation  21  with  E of  equation  20  yields 

-n  (7Q)  x 7'  (7q)  = n (7q)  x JJ  jwp0  7s  (7)  . 7 (7,  7q)  dA 

(Eq.  22) 

For  the  thin-wire  approximation,  limiting  attention  to 
circular  cross-section  bodies  of  diameters  small  compared  with  the 
wavelength,  the  azimuthal  current  may  be  neglected,  and  equation  22 
becomes 
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-n  (rQ)  x E1  (r( 


' n (ro>  * if  J J { K Js  <r)  [ * + ] 9 <r.  ro)  } 


(Eq.  23) 

where  s is  the  unit  tangent  vector  at  r pointing  in  the  direction  of  the 

current.  A scalar  integral  equation  for  the  current  is  obtained  by 

taking  the  dot  product  of  equation  23  with  the  unit  tangent  vector  sq 

at  tie  observation  point  r as 

o 


‘V^^O*  = ^ \\  jU,M0  JS  (ro)  [s-s0+  (s* 


v)  (s  . v)  — , g (r,  r ) dA 
0 , / o 


(Eq.  24) 


If  the  assumption  is  now  made  that  is  independent  of 
the  azimuthal  variable,  equation  24  can  be  written 


rl  . - 1 

-s  • E (r  ) = t — 
o o 4tt 


j awu  J$  (s) 


r 7 
f I*  1 3 

J |_S  • So  ‘ ,2  3s3so 


g(r,  r ) d<t»ds 


(Eq.  25) 


where  a is  the  wire  radius  and  the  s integration  is  over  the  entire 
length  of  wire  L.  A final  approximation  is  that  the  current  may  be 
realistically  represented  as  a filament  of  strength  1^  (s)  = 2Tta  (s) 
flowing  on  the  wire  axis  while  the  field  is  evaluated  on  the  wire 
surface,  allowing  equation  25  to  be  written 

-s  • E 1 (7o>  = (‘ja)MoAir)  \ 's  (s)  [S‘  so  ' 7T  373T-]9(7*  ?o)  ds 

O " L k 0-1 

(Eq.  26) 


I 


t 


t 


where  |r  - rQ|  is  now  measured  from  the  wire  axis,  or  source  point,  to 
the  observation  point  on  the  surface,  which  can  thus  never  be  closer 
than  the  wire  radius  a.  By  considering  the  current  as  a tubular  sheet 
centered  on  the  wire  axis  while  evaluating  the  electric  field  at  the 
wire  axis,  one  can  resolve  the  ambiguity  in  the  azimuth  involved.  The 
form  of  equation  26  is  not  changed  using  this  convention,  but  the 
interpretation  of  the  tangential  field  evaluation  is  simplified  when 
nonparallel,  nonplanar  wires  are  considered. 

The  thin-wire  approximation  which  leads  from  the  electric 
field  integral  for  a surface  current  distribution  to  equation  26 
involves  the  assumption  that  the  wire  radius  a <<  X so  that:  (l)  azimu- 

thal current  flow  around  the  wire  may  be  neglected;  (2)  the  longitudinal 
current  is  independent  of  azimuth  and  may  be  represented  as  a filament 
along  the  wire  axis;  and  (3)  that  the  surface  integration  can  be  replaced 
by  a line  integration  along  the  wire  at  a radial  distance  a^  so  that  the 
minimum  source-to-observation  point  distance  is  an  thus  avoiding  the 
singularity  in  the  kernel  of  the  integral  which  would  occur  at  r = r . 

This  thin-wire  approximation  has  been  applied  to  radiation 
and  scattering  problems  with  a great  deal  of  success.  Until  fairly 
recently,  for  example,  linear  antenna  theory  was  almost  exclusively 
restricted  to  the  thin-wire  approach;  the  same  observation  applied  to 
scattering  from  finite  length  cylinders.  At  the  same  time,  the  antenna 
and  scattering  solutions  were  largely  confined  to  dipoles  on  the  order 
of  a few  wavelengths  long.  This  results,  for  the  most  part,  from  the 
approximate  analytic  approaches  required  due  to  the  difficulty  involved 
in  obtaining  a numerical  solution  of  the  integral  equation.  While  Wu 
(I960)  and  Chen  ( 1 968 ) extended  the  analytic  methods  developed  by  Hallen 
(1938)  and  King  (1956)  to  antenna  and  scatterers  without  this  length 
restriction,  their  results  are  quite  complicated  in  form  and  in  addition 
have  not  been  subjected  to  extensive  experimental  comparison. 
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The  reason  for  concentration  on  analytic  solutions  to  the 
integral  equation,  until  the  past  few  years,  was  the  lack  of  sufficiently 
powerful  computers  to  provide  the  capability  to  obtain  a completely 
numerical  solution.  With  the  present  development  of  both  high-speed 
computers  and  advanced  methods  of  numerical  analysis  this  is  no  longer 
the  case:  significant  progress  has  been  made  in  extending  accurate 

numerical  solutions  to  more  complex  geometries.  In  the  following  section, 
the  method  of  solution  will  be  outlined  for  the  electromagnetic  properties 
of  structures  to  which  the  thin-wire  integral  equation  is  applicable, 
b.  Reduction  to  a Linear  System  (Collocation) 

A numerical  solution  to  an  integral  equation  may  perhaps 
be  best  undertaken  using  the  method  of  moments.  This  is  a well-founded 
mathematical  technique  for  finding  the  unknown  by  forcing  the  integral 
equation  to  be  satisfied  in  some  prescribed  fashion  over  the  range  of 
the  integral  operator.  GEMACS  is  based  on  the  thin-wire  electric  field 
integral  equation. 

Equation  26  may  be  written  symbolically  as: 


£(f)  = g 


(Eq.  27) 


following  Harrington's  (1968)  notation.  The  solution  of  equation  26 
(or  of  equation  27)  is  obtained  by  the  method  of  moments.  An  intuitive 
approach  to  solving  equation  27  for  the  unknown  function  f i s to  set  f 
equal  to  a constant  f.  within  N subintervals  of  the  domain  of  £ , and  to 
require  equation  27  to  be  satisfied  at  N points  over  the  range  of  £, 
obtaining  N equations  in  the  f.  unknowns.  This  is  a specialized  applica- 
tion of  the  method  of  moments  which  is  more  generally  written  as  follows. 
Let 

f = ^ a f 

with  the  basis  function  f defined  in  the  domain  of  £ so  that  equation 

n 

27  may  be  wr i tten 
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Lan  ^ (fn)  = 9 


(Eq.  28) 


Then,  with  the  set  of  weighting  functions  w^,  defined  in  the  range  of 
£ , the  inner  product  denoted  by  < > is  formed  as 


^an  < V ^ (fn} 


V 9 


(Eq.  29) 


where  m = 1,2,3  ....  Equation  29  can  be  written  in  matrix  form  as 


[G  ] [a]  = [S  ] 
mn  n m 


where 


G =<  wm.  £ (f) 

mn  \ m n 


S = / w , g 
m \ m 3 


and  the  matrix  G is  referred  to  as  the  structure  matrix.  If  the 
mn 

inverse  of  G exists,  then  the  a can  be  found, and  thus,  the  function  f, 
mn  n 

which  is  the  desired  solution,  for  any  specified  source  function  S . 

m 

The  proper  choice  of  weight  functions  and  basis  functions, 
as  well  as  the  subsectioning  of  the  domain  of  £ , is  not  an  obvious  one. 
Although  there  is  some  leeway  in  the  matter,  careful  consideration  of 
the  physics  of  the  problem  and  the  nature  of  the  expected  solution  will 
show  that  some  representations  for  the  f will  be  more  efficient  than 
others  in  terms  of  computer  time  and  accuracy.  Constant,  linear,  quadra- 
tic, tr igonometr ic, and  Fourier  series  have  all  been  used  for  this  role. 
The  weighting  functions  have  generally  been  more  restricted  in  choice 

than  f . The  special  case  w = f is  referred  to  as  Galerkin's  method, 
n n n 


More  often,  the  weights  are  6-functions,  a method  referred  to  as 

collocation,  so  that  the  inner  product  equation  25  merely  becomes 

the  sequence  of  values  Z (f  ) and  g . These  are,  respect i ve 1 . , 

n m m 

the  tangential  electric  fields  due  to  current  segment  n at  observa- 
tion point  m and  the  tangential  incident  electric  field  at  observa- 
tion point  m. 

c.  Current  Expansion 

The  GEMACS  program  employs  the  collocation  method  with 
constant,  sine  and  cosine  terms  for  the  f^  segment  or  current  function, 
i .e. , 

N N 

I (s)  = ^ U.(s)[A.  + B.  sin  k(s-s.)  + C.cos  k (s-s.)]  = U . (s)  I . (s) 

^ J J J J J J ^ J J 

j = 1 j - 1 

(Eq.  31) 

where  U.(s)  is  1 when  s i s on  segment  j and  zero  otherwise.  Equation 
31  appears  disadvantageous  because  three  constants  are  required  to 
specify  the  current  on  each  segment,  so  that  apparently  3N  linear  equa- 
tions need  be  solved;  however,  it  is  not  necessary  to  employ  the  integral 
equation  itself  to  find  the  extra  unknowns  introduced  by  the  sinusoidal 
expansion.  Two  of  the  three  constants  for  each  segment  may  be  obtained 
by  requiring  the  current  in  adjacent  segments  to  satisfy  some  specified 
mutual  conditions.  In  the  GEMACS  thin-wire  program,  the  extrapolated 
current  from  a given  segment  is  forced  to  match  the  center  current 
values  in  two  adjacent  segments  to  satisfy  the  required  condition  for 
two-wire  junctions  as  follows. 

Let  the  current  on  segment  j be  expressed  as 

I . (s)  = A.  + B.  sin  k(s-s.)  + C.  cos  k (s-s.) 

J J J J J J 

(Eq.  32) 
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with  Sj  the  midpoint  coordinate  (Yeh  and  Mei,  1967).  Also,  let  segment 

j be  connected  to  segments  j-1  and  j+1  at  its  minus  and  plus  reference 

ends  respectively  with  the  reference  directions  on  all  three  segments 

the  same.  Evaluation  of  I . at  s.,  s.  ,,  and  s.  . results  in 

J J J-1  J+1 


A.  + C.  = I . 
J J J 


A . + B . s . . . + C . c . = I . 

J J J-l.J  J J “ 1 » J J' 


(Eg.  33) 


A.+B.  s.  , . + C . c . , . — I 

j j j+i,j  j j+i,j  j+i 


where  dj  is  the  length  of  the  j segment  and 


k (dj  + 1 + d 


j,/2] 


c j+1 , j 


Solution  for  A.,  B.,  and  C.  in  terms  of  I.  I.,  and  I.  , provides  an 
J J J J-1  J J+1 

equation  of  the  form 


Ij(s')  - x.U')!..,  ♦ V.(s')  i.  * z.  (s')  ij+, 


(Eq.  34) 


where  X.,  Y.,  Z.  contain  the  coefficients  A.,  B.,  C..  The  system  of 
J J J J J J 

equations  which  results  from  the  collocation  solution  to  the  integral 

equation  is  thus  seen  to  involve  as  unknowns  the  N current  samples  at 

the  centers  of  the  N segments  into  which  the  structure  is  divided. 

✓The  extension  of  the  interpolation  procedure  to  multiple 

junctions  is  straightforward.  Consider  the  case  where  segment  j is 

connected  to  m segments  numbered  j + 1,  , j + m at  its  plus  end  and 
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the  single  segment  j-1  at  its  minus  end.  Then  only  the  equation  repre- 
senting equation  32  evaluated  at  sj+j  is  modified  and  becomes 


A.  + — V'  (B.  s.  . + C.  c.  .)  = — V'  I.  (Eq.  35) 

J m LJ  J i »J  J i , J m ' j v H 

i=j+l  i=j+l 

which  comes  from  interpolating  I . to  the  midpoints  of  the  m connected 

wires  where  it  is  equated  to  the  average  midpoint  value.  A solution  for 

the  A.,  B.  and  C.  in  terms  of  the  midpoint  currents  I.  ,,  I.  ,,  .... 

J J J J-1  j + r 

I.  follows  as  before.  A multiple  junction  at  the  minus  end  of  the 
j+m  J 

segment  is  similarly  treated. 

The  sinsusoidal  current  expansion  appears  to  make  the 
system  of  equations  resulting  from  collocation  somewhat  more  involved. 
But  the  required  computer  time  is  not  significantly  increased  when 
compared  with  the  same  number  of  current  unknowns  without  using  the 
sinusoidal  expansion.  Other  current  expansion  funct ions- 1 i near , quadra- 
tic, Fourier  series-could  be  used  in  place  oF  the  constant-sine-cosine 
expression,  but  this  particular  expansion  has  a number  of  additional 
advantages  over  the  other  possibilities  mentioned.  For  instance,  a 
solution  for  the  current  to  a specified  accuracy  for  a half-wave  dipole 
scatterer  and  antenna  requires  the  fewest  current  segments  using  the 
sinusoidal  expansion  (Neutreuther , et  al.,  1968).  This  advantage  would 
be  expected  to  carry  over  to  more  complex  geometries.  Second,  the 
solution  will  more  accurately  exhibit  the  required  dependence  on  wire 
radius  (Andreason,  1968)  because  the  constant  current  term  produces 
infinite  tangential  electric  field  on  the  current  axis,  as  opposed  to 
the  sine  and  cosine  terms  which  do  not. 

Third,  the  parallel  and  perpendicular  electric  field 
components  (due  to  the  sine  and  cosine  current  terms)  and  the  tangential 
field  components  (due  to  constant  current  terms)  may  be  analytically 
evaluated.  This  eliminates  the  necessity  for  extensive  numerical 
integration  to  evaluate  all  the  elements  in  the  coefficient  matrix 
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G . Only  the  perpendicular  electric  field  excited  by  the  constant 
current  terms  requires  numerical  integration  and  this  is  handled  by 
applying  the  Romberg  variable- interval  width  technique  to  the  difference 
integrand.  This  is  discussed  in  the  appendix  of  the  AMP  Engineering 
Manual . 

d . Calculation  of  the  Structure  Matrix 

The  form  of  the  matrix  elements  which  result  from  applying 
the  method  of  collocation  to  equation  26  is  considered  in  the  following 
discussion.  Each  entry  G.^  in  the  structure  matrix  represents  the 
tangential  electric  field  at  observation  point  i on  the  structure  produced 
by  unit  current  flowing  on  segment  j.  The  boundary  condition  on  the 
tangential  electric  field  is  enforced  at  each  observation  point.  The 
collocation  method  of  solving  the  integral  equation  is  thus  basically 
one  of  calculating  electric  field  components  at  specific  points  due  to 
the  current  induced  on  the  structure. 

The  thin-wire  approximation  involves  the  explicit  assump- 
tion that  the  effects  of  azimuthal  currents  can  be  neglected  in  compari- 
son with  those  of  axially  directed  currents  and  that,  in  addition,  the 
cylindrical  tube  of  axial  current  has  no  azimuthal  dependence.  The 
former  assumption  allows  one  to  consider  only  one  current  component 
rather  than  two,  while  the  latter  provides  partial  justification  for 
reducing  the  surface  integral  to  a line  integral.  It  may  be  deduced 
from  an  examination  of  equations  25  and  26,  however,  that  even  where 
is  independent  of  if^the  kernel  of  the  integration  equation  depends  in 
general  upon  both  4>  and  s.  However,  the  integrand  is  independent  of  <p 
in  the  special  case  where  the  observation  point  is  located  on  the  axis 
of  a linear  tube  of  current,  and  the  <j>  integration  of  equation  25  may 
be  replaced  without  approximation  by  the  factor  2ir. 

Consequently,  the  observation  points  are  located  where 
the  tangential  electric  field  is  to  be  calculated,  on  the  axis  rather 
than  on  the  surface  of  each  wire  segment.  The  <f>  integration  in  equation 
25  is  thus  exact  for  the  self-field  as  well  as  the  mutual  fields  for 
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all  current  segments  having  a common  axis.  In  addition,  the  possible 
ambiguity  involved  in  evaluating  the  incident  field  over  a 2tt  variation 
in  <J>  on  the  wire  surface  is  resolved.  As  a final  point,  the  observation 
point  is  always  at  least  as  far  as  the  wire  radius  from  the  source 
point. 

When  the  mutual  fields  of  nonaxial ly  aligned  current 
segments  are  required,  the  <p  integration  is  not  so  simply  performed. 

And,  if  no  approximation  were  used,  the  <J>  integration  would  require 
numerical  evaluation.  The  most  obvious  approach  is  to  consider  the 
tubular  current  source  to  approximate  a linear  filament  on  the  wire 
axis,  a procedure  which  again  replaces  the  <J>  integration  by  a 2tt  factor. 
Unfortunately,  this  approximation  eliminates  the  influence  of  the  wire 
radius  from  all  mutual  field  terms  on  the  phase  change  and  geometrical 
attenuation  of  the  field  caused  by  the  separation  of  the  source  and 
observation  points. 

An  alternative  to  the  above  method  is  replacement  of  the 
current  tube  by  a current  filament  which  is  not  located  on  the  wire  axis 
but  is  displaced  in  distance  from  it  by  the  wire  radius.  The  direction 
of  displacement  is  perpendicular  to  the  plane  of  the  wire  axis  and  the 
line  joining  the  observation  point  and  wire  axis  midpoint  (the  observa- 
tion point  for  the  self  term  field).  The  geometry  of  this  method  is 
shown  in  figure  1.  The  improved  agreement  obtained  between  theory  and 
experiment  for  a 16-sided  regular  polygon  backscatter  cross  section 
demonstrates  that  this  approach  is  more  realistic. 

To  summarize  briefly,  the  surface  integral  is  reduced  to 
a line  integral  by  neglecting  azimuthal1  currents  and  azimuthal  variation 
of  the  axial  current.  Self-field  terms  are  calculated  with  the  observa- 
tion point  on  the  axis  of  a cylindrical  current  tube,  while  mutual  field 
terms  are  calculated  at  the  same  observation  point  with  the  current 
represented  as  a filament  displaced  from  the  wire  axis  by  the  wire 
rad i us . 


^aken  here  to  mean  the  direction  measured  along  the  intersection  of 
the  current  tube  surface  with  a plane  perpendicular  to  the  axis  of 
the  current  tube. 
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OBSERVATION  SEGMENT 


1 


i 


(b ) Geometric  parameters  for  field  evaluation 


Figure  1.  Geometry  of  Current  Filament 
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form 


The  integral  equation  26  can  now  be  written  in  the 


N 

- a™  /mv  c -»L- 

o J m n ^2  3s3s 

n=  1 *" 


g(r,  r)  I (s)  ds 

m n 


(Eq.  36) 
m = 1 , . . . , N 


where  AS^  denotes  the  length  of  segment  n,  is  the  tangential  compo- 
nent of  the  incident  electric  field  at  the  mth  segment,  and 


9(rm’  r)  = 


-jk  \Jp2  + 


( s - s)‘ 
m 


I 2 A , v 2 

P + (s  - s) 
m 


(Eq.  37) 


It  should  be  noted  that  the  integration  over  L has  been  reduced  to  a 
summation  of  N separate  straight-wire  segment  integrals.  It  is  conve- 
nient to  rewrite  equation  36  ;n  terms  of  cylindrical  coordinates 
referred  to  the  wire  segment  being  integrated.  Then  one  obtains 


E = (jwu  Mir) 
m o 


N 

2[  * 1 a 

3 Zn  ■ Sm  2 Sz  3s 

[_  k n m 


7 r — g(r  ,z)l  (z  ) dz 

oz  3s  mn  n n n n 

n m 


where 


(Eq.  33) 


g(r  , z ) = exp  (- i kr  )/r  (Eq.  39) 

mn  n r mn  mn  ^ 


r = \/(z  - z )2  + p 2 + a2 

mn  \/  mn  n mn  n 


(Eq.  1,0) 


are  the  radial  and 


a is  the  radius  of  wire  segment  n,  and  p and  z 
n mn  mn 

z-coord i nates  of  the  observation  point  at  the  center  of  segment  m referred 
to  the  midpoint  of  segment  n,  as  shown  in  figure  1. 

e.  Impedance  Loading 

The  discussion  has  thus  far  been  limited  to  the  case  of  a 
perfectly  conducting  scatterer.  The  approach  may  be  generalized  to 
allow  for  loading  of  the  structure  by  introducing  a voltage  drop  term  in 
the  integral  equation.  If  the  impedance  loading  per  unit  length  on 
segment  m is  Z^,  then  equation  38  becomes 

E * - I Z = same  right-hand  side  ( Eq . 41) 

m mm 

f . Current  Solution 

Having  evaluated  the  mutual  impedance  elements  (the 
structure  matrix),  equation  26  can  be  written  in  matrix  notation  as 

N 

y]  Gjj  fj  - -e!  i = 1,  2,  ...  N (Eq.  42) 

J-l 


where  G is  the  structure  matrix,  I.  is  the  unknown  current  at  wire 

I J 

segment  j,  and  E.  is  the  incident  tangential  field  at  segment  i. 
Equation  42  is  solved  in  the  form 


N 


I . = 
J 


E. 


(Eq.  43) 


J-l 

The  operation  implied  by  this  equation  may  be  accomplished  via  inversion, 
factorization  of  the  G matrix,  or  by  iteration.  In  the  collocation 
solution  of  the  wire-antenna  problem  used  in  the  GEMACS  program  for  arbitrary 
thin-wire  s tructures, the  solution  of  equation  43  for  the  current  will 
represent  the  major  portion  of  computer  program  execution  time  for 
complex  structures  containing  a large  number  of  segments.  It  is,  therefore, 
of  the  utmost  importance  to  have  efficient  solution  procedures  available 
for  this  tvDe  of  structure. 
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. Solution  Process  Implementation 


The  basic  solution  process  involves  decomposing  the  structure 
matrix  [G]  (or  that  portion  of  it  extracted  for  a BMI  solution),  per- 
forming forward  elimination  and  back  substitution,  and  performing  matrix 
multiplication.  The  performance  of  these  functions  is  straightforward 
when  „ H of  the  data  reside  in  core.  However,  GEMACS  was  designed  to 
handi-=‘  large  problems  and  thus  core  resident  routines  are  not  very 
efficient.  The  primary  function  is  that  of  decomposing  a matrix  into  an 
upper  and  lower  triangular  matrix  when  none  of  the  matrices  will  fit  in 
core  storaqe.  The  algorithm  employed  in  GEMACS  is  specifically  tailored  to 
the  data  storage  method  for  matrices.  All  matrices  are  stored  by  column, 
that  is,  each  column  requires  a separate  read  and  write.  For  this 
discussion,  a column  will  be  referred  to  as  a record.  While  this  is 
obviously  not  the  best  1/0  scheme,  it  is  likely  to  be  fairly  good  since 
larger  problems  typically  require  several  hundred  elements  and  each 
record  will  contain  twice  that  number  of  entries  since  the  matrix  is 
complex.  Therefore,  the  system  buffer  for  I/O  will  be  quite  adequately 
used  with  a single  record.  Enlarging  the  buffers  usually  increases  the 
core  size  and  will  not  significantly  reduce  run  time.  In  addition,  all 
I/O  is  in  ANSI  FORTRAN  IV  for  compatibility  with  a large  number  of 


mach i nes . 


The  function  of  the  solution  processor  is  to  find  the  solution 


I to  the  set  of  simultaneous  linear  equations: 

G I = V 


(Eq . 4A) 


The  method  employed  is  known  as  lower/upper  triangular  decomposition 
which  is  a variant  of  Gaussian  elimination.  The  matrix  G is  decomposed 
into  lower  (L)  and  upper  (U)  triangular  matrices  such  that 


G = L U 


G I = V 
LUI  = V 

Ul  = L_1  V = Y 
I = u” 1 Y 

* Note  that  the  inverses  of  L and  U are  never  found  and  the  notation  is 

used  to  show  the  solution  logic  only.  Also,  the  inverse  of  G is  not 
found  unless  explicitly  requested. 

When  decomposing  a matrix,  one  proceeds  down  the  diagonal 
modifying  all  elements  below  and  to  the  right  of  the  diagonal  by  opera- 
tions performed  upon  elements  of  the  diagonal  row  and  column.  Thus, 
once  the  i**1  diagonal  element  has  been  used,  all  column  elements  below 
and  row  elements  to  the  right  of  the  diagonal  will  not  be  referenced 
, again.  These  are  the  elements  of  the  lower  and  upper  triangular  matrices 

and  may  be  written  out  immediately  to  their  respective  files.  Likewise, 
the  elements  of  the  square  submatrix  remaining  may  be  written  to  a 
peripheral  file  and  when  all  of  the  elements  have  been  processed,  all 
» future  elements  reside  in  this  square  submatrix  which  is  smaller  by  one 

row  and  column  than  the  source.  This  procedure  is  repeated  N- 1 times 
where  N is  the  dimension  of  the  original  matrix.  In  this  way,  the  point 
is  reached  where  the  entire  remaining  submatrix  will  fit  in  core  and  may 
be  decomposed  using  normal  codes.  This  procedure  is  illustrated  in 
figure  2 for  the  first  three  diagonal  elements  of  a 10  x 10  matrix. 

Note  that  even  though  the  elements  of  the  upper  triangular 
matrix  are  row  elements,  they  will  be  stored  as  column  elements.  This 
will  be  compensated  for  during  the  elimination  process. 

The  types  of  decomposition  encountered  are  referred  to  as  row 
or  column  decomposition  depending  on  the  sequence  of  operations  performed. 
GEMACS  employs  column  decomposition  and  diagonal  pivoting;  that  is,  no  row 
or  column  interchanges  take  place.  The  lack  of  interchanging  rows  or 


Zh 


N = 10 


n = 1 


COLUMN  ) 
OF  LOWER 


N = 9 


COLUMN  2 
OF  LOWER 


N = 8 


•COLUMN  I OF  UPPER 


9X9  SUBMATRIX  FOR  NEXT  ROUND 


’COLUMN  2 OF  UPPER 


'3X3  SUBMATRIX  FOR  NEXT  ROUND 


COLUMN  3 OF  UPPER 

7X7  SUBMATRIX  FOR  NEXT  ROUND 


Figure  2.  Illustration  of  Matrix  Reduction  During  Decomposition 
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columns  may  lead  to  problems  for  ill-conditioned  matrices;  however,  the 
matrices  generated  by  GEMACS  do  not  display  this  characteristic  except 
possibly  near  high-Q  resonances  which  are  unikely  on  large  structures. 

The  matrix  decomposition  is  monitored  to  detect  the  occurrence  of 
instability  or  errors  due  to  round-off  accumulation.  Such  indication 
did  not  occur  during  the  extensive  development  work  on  the  BMI  scheme 
for  a wide  variety  of  shapes.  The  incorporation  of  a pivoting  capability 
in  GEMACS  is  possible  if  a need  is  demonstrated. 

The  algorithm  for  GEMACS  decomposition  at  the  round  is: 

G ( J , I ) = G (J , 1 )/G (J , J ) for  I = J+l,  N (Eq-  ^5) 

G (K,  I ) = G(K,  I ) - G(J,I)/VG^,J)  for  K = J+l,  N (Eq.  A6) 

Note  that  the  second  term  in  equation  A6  contains  the  con- 
stant G (J , I )/G  (j , J ) for  I and  J fixed.  Then  as  K varies,  contiguous 
elements  of  the  I ^ column  are  modified.  In  the  FORTRAN  code,  this 
permits  simple  subscripts  to  be  used  with  incremental  steps  and  alleviates 
the  need  for  determining  the  storage  address  of  the  elements  for  each 
value  of  I.  This  results  in  a more  efficient  code  which  executes  much 
faster  than  codes  which  perform  decomposition  by  rows.  Also,  this 
method  makes  it  computationally  simpler  to  decompose  real  or  complex, 
banded  or  nonbanded  matrices  using  the  same  subroutine. 

The  result  of  the  decomposition  is  lower  and  upper  triangular 
matrices  written  by  record  to  two  separate  peripheral  files.  The  matrices 
could  have  been  combined,  however,  they  are  never  both  needed  simultane- 
ously. Therefore,  more  of  each  will  fit  into  the  available  core  using 
this  method.  Normally  one  of  the  matrices  will  have  ones  (1)  on  the 
diagonal  while  the  other  will  have  the  diagonal  elements  as  modified 
by  decomposition.  The  GEMACS  code  places  the  diaqonat  elements  on  both 
matrices  since,  if  the  original  matrix  was  transposed  (as  in  the  case 
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for  the  GEMACS  structure  matrix)  the  role  of  the  lower  and  upper  matrices 
is  interchanged.  If  L and  U are  the  lower  and  upper  matrices  resulting 
from  decomposition  of  G,  then 


I 

i 

i 


G = LU 


gt  = (lu)t 


GT  = LT  UT 


(Eg.  47) 


(Eq.  48) 


where  LT  and  UT  are  the  lower  and  upper  triangular  matrices  obtained 
from  transposing  U and  L.  Therefore,  the  role  has  reversed  and  the 
logic  has  been  incorporated  in  the  GEMACS  code  to  always  assume  the  lower 
triangular  matrix  has  unit  diagonal  elements. 

Using  the  algorithm  described  above,  GEMACS  needs  room  for  only 
three  columns  of  the  matrix  in  core  in  order  to  perform  the  decomposition 
(2  columns  for  decompositon  and  1 column  to  accumulate  the  elements  of  the 
upper  triangular  matrix).  The  lower  matrix  is  stored  on  a scratch  file. 

Once  the  matrix  has  been  decomposed,  the  data  to  be  used  in 
the  lower  matrix  are  recovered  by  column  until  the  available  core  storaqe 
is  used.  Forward  elimination  is  the  process  of  solving  the  system 

LY  = V (Eq.  49) 


where  V is  the  original  right-hand  side  and  Y is  originally  the  null 
vector.  Storage  for  Y is  not  actually  required  since  V may  be  overwritten. 

The  elements  of  V are  modified  according  to 

i-1 


j=i 


l V. 

U J 


(Eq.  50) 
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where  is  an  element  of  the  lower  triangular  matrix.  It  is  seen  that 
the  element  v.  depends  only  on  the  elements  v^.  which  precede  it,  thus  the 
name,  forward  elimination.  Therefore,  v.  may  be  accumulated  for  those 
rows  of  V for  which  the  columns  of  L are  in  core.  Note  that  the  elements 
of  L were  written  out  as  a column  even  if  G is  a transposed  matrix  since 
the  rows  of  U become  the  columns  of  L. 

Once  Y is  found,  the  system  that  remains  is 


U I = Y 


(Eq.  51) 


where  U is  the  upper  triangular  matrix  with  1 row  per  record,  I is  the 
solution  vector,  and  Y is  L 'v  or  a modified  right-hand  side.  Since  U is 
a lower  triangular  matrix,  one  starts  at  the  bottom  and  works  back  up 
the  right  hand  side.  The  elements  of  I again  overwrite  Y or  V and  are 
given  by 

N 

vi  ■ (vi  • S uij  vj)/uii  (E<-  52) 

j=i+l 


where  u.^  is  an  element  of  the  upper  triangular  matrix.  Since  the  last 
elements  of  U are  needed  first,  the  GEMACS  code  will  determine  how  many 
rows  will  fit  and  fetch  the  data  in  the  proper  sequence.  Since  ANSI 
FORTRAN  does  not  support  random  access  I/O,  this  can  be  a time  consuming 
and  expensive  process  for  very  large  matrices.  This  expense,  as  well  as 
the  expense  of  decomposition,  is  considerably  reduced  for  banded  matrices, 
however,  when  using  BM  I , the  solution  process  is  repeated  for  each  itera- 
tion. A detailed  discussion  of  BM I is  presented  in  section  D. 

Implicit  in  using  BM I is  the  matrix  multiplication  involved  in 
finding  the  RHS  (right-hand  side)  at  each  iteration,  i,  prior  to  the 
solution  of  equation  A9. 


RHS.  = V - (GL  + GU)  I M 


(Eq.  53) 
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The  elements  of  equation  53  are  matrices  where  GL  and  G11  are  those  parts 
of  the  structure  matrix  not  included  in  the  band.  V is  the  excitation 
and  I | is  the  last  solution  obtained.  Remembering  that  we  have  core 
storage  for  at  least  three  columns  of  the  structure  matrix,  V.,  |. 
and  one  column  of  G wi 1 1 fit  in  core.  The  elements  of  GL  and  GU  remain 
imbedded  in  G.  The  elements  of  RHS  for  a bandwidth  of  M are  (where  i 
and  j are  subscripts  of  the  matrices): 

N 

RHS  =Vi-  E (GL  + GU  ) I ' (E  53(a)) 
j < i -m  J J 1 M 

j > i+m 

and  are  accumulated  as  the  blocks  G are  read  into  core.  The  same  tech- 
nique of  partial  sum  accumulation  may  be  used  to  multiply  matrices 
together  which  will  not  fit  into  core,  and  it  is  even  more  efficient  than 
conventional  methods  when  they  do  fit, since  again,  the  innermost  FORTRAN 
DO  loops  reference  continuous  data, and  thus, an  optimizing  compiler  will 
not  generate  the  indexing  locations  code. 

In  discussing  how  to  handle  large  matrices,  the  idea  of  very 
large  core  machines  usually  comes  up.  While  such  machines  exist,  few 
can  directly  access  very  large  arrays.  For  instance,  while  the  CDC  7600 
may  have  1 million  words  of  l.CM  (Large  Core  Memory),  only  131,000  words 
may  be  addressed  in  any  single  array.  This  is  due  to  the  size  of  the 
address  word  used  by  the  compiler  in  indexing  computations.  It  is 
highly  unlikely  that  anyone  is  going  to  be  willing  to  dedicate  more 
than  1 8 bit  words  to  address  storage, and  thus, 131,000  is  the  limit. 

Without  assembly  language  routines,  this  necessitates  the  use  of  out-of- 
core techniques  and  GEI1ACS  has  been  specifically  designed  for  this  within 
the  constraints  imposed  by  ANSI  FORTRAN  IV. 

5 . Observables  Computation 

The  GEMACS  code  will  compute  the  near  and  far  electric  field 
upon  request.  These  quantities  are  computed  exactly  as  in  the  AMP  code. 

Portions  of  the  AMP  Systems  Manual  describing  the  AMP  subroutines  EFLD 
and  FFLD  are  reproduced  here.  The  GEMACS  code  does  not  include  the  ground 
wave  capability  present  in  the  AMP  code. 

fl 
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The  far  electric  field  due  to  line  currents  can  be  written 


E (r)  = jwu 


[ (k (k  * jeJk*r  T(?,)  djl))  • j eJk* 


I (r 1 ) d* 

(Eq.  5*0 

where  r is  the  position  vector  of  the  observation  point,  r'  is  the  posi- 
tion vector  of  the  source  point,  k is  in  the  direction  of  propagation 
with  a magnitude  of  2tt/A.  Specialized  to  straight  wire  segments  as  used 
in  the  GEMACS  formulation 

T ,-s  . \ e"jkr  V jk*  R.  T 

E <r)  ■ 1 v T7 r ,L  e 1 L 


k (k  • Q. ) 


2tt 


i = l 


(Eq.  55) 


where  n is  free  space  impedance,  R.  is  the  position  vector  of  the 

^ th  ^ 

center  of  the  i segment  and 


S/2 


-“i  $ 

-(S/2) 


e j 2tt ( k • u.)t  I . (t) 


dt 


where  u.  = cos  a.  cos  g.  x + cos  a.  sin  (5.  y + sin  a.  z which  is  the 

i i i t h i i i 

reference  direction  of  the  i segment  with  the  angles  defined  as  shown 

below,  tu.  = r'/A  -R./A,  and  s is  the  segment  length  in  wavelengths. 
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— ----- — 


Wi  th 


I . (t)/X  = A.  + B.  sin  2vt  + C.  cos  2irt 
i ii  i 


integration  of  Q.  yields 


„ f sin  it 

= u.  A.  

I |_  I TTW 

(sin  it  ( 1 

177TT 


TTW.  S 

I 


in  7T ( 1 + w.)s  sin  ir(l  - w.)s 


+ J 


(sin  tt  ( 1 + 
2tt  ( 1 + w . 


(sin  it  ( 1 + w.)s  sin  tt  ( 1 - w.)sl 
2tt  ( 1 + w. ) + 2tt  ( 1 - w . ) ) 


where  w.  = - k • u 


2u ( 1 - w. , 


(Eq.  56) 


(Eq.  57) 


Note,  the  term  k (k  »Q.)  in  equation  55  is  completely  radial  and 
cancels  the  radial  component  of  Q. . This  term  is  ignored  in  GEMACS  since 
the  desired  transverse  components  will  be  computed  by  a dot  product. 

Thus  for  program  use  only  and  with  the  understanding  that  only  transverse 
components  will  be  used,  we  write 


n -jkr  r — i 

E(r)  = -j  ^ hjr  2 


jk  • R. 


(Eq.  58) 


Ground  effects  are  included  by  means  of  an  image  and  the 
appropriate  reflection  coefficients.  The  z component  of  the  segment; 
reference  direction  vector  u changes  sign  for  the  image  as  shown 
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and  0 is  measured  from  z. 

For  a semi - i nf i n i te  ground,  the  electric  field  of  the  image 
can  be  calculated  from  equation  55  and  substituted  into  equation  59 
to  calculate  the  reflected  field. 

To  compute  the  near  field  at  a point  in  space  due  to  the 

current  on  a wire  structure,  the  field  is  computed  for  three  current 
distributions:  sine,  cosine,  and  constant  functions  on  all  wire  seg- 

ments and  then  summed  vectorially. 
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The  wire  segment  is  considered  to  be  located  at  the  origin 
of  a local  cylindrical  coordinate  system  with  the  point  at  which  the 
field  is  computed  being  (p1,  <J>',  z').  The  geometry  for  a filament  of 
current  of  length  A is  shown  below. 


Segment 


For  a sine  or  cosine  current  d i str i but  ion, the  field  can  be  written  in 


1 

closed  form.1 

The  p arid  z 

field 

components  for  a 
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Stratton,  J.  A.,  Electromagnetic  Theory,  McGraw  Hill  Book  Co.,  New  York, 

1941,  p.^5^. 
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A/2 


Ez  (p'>  Z,)  = T 


(-j  -f-  ) 


e‘jkr2 


(kr2)‘ 
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-jkr 


-A/2 


-0  + ) (kz1  + kA/2)  e 


-jkr 


(kr1 )' 


dz  + (j  + — — ) (kz1  -kA/2) 


(Eq.  63) 


These  expressions  are  separated  into  real  and  imaginary  parts  for 
evaluation  in  the  program.  The  coordinate  p'  for  a wire  segment  is  taken 
as  the  distance  from  the  observation  point  to  a point  on  the  side  of  the 
segment  as  shown  below. 


Also,  the  components  E are  multiplied  by  p/p'  to  account  for  the  change 
in  vector  direction. 

Ground  image  contributions  are  taken  into  account  in  the  same 
manner  as  for  the  far  field  computation. 
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BANDED  MATRIX  ITERATION  THEORY  AND  DEVELOPMENT  SUMMARY 


1 


. 1 . Introduction 

This  section  is  a summary  of  the  theory  and  development  of  the 
j BMI  (Banded  Matrix  Iteration)  solution  technique  for  the  linear  simulta- 

neous equations  arising  from  the  thin-wire  method  of  moments  formalism. 
t Standard  methods  of  computer  solution  are  too  expensive  for  application 

to  large  problems.  The  BMI  technique  was  developed  to  reduce  this  cost. 

It  was  chosen  from  a variety  of  possible  new  alternative  methods  after  a 
review  of  the  literature  (reference  1).  The  first  application  was  for 
* single  straight  wires,  either  antennas  or  scatterers,  up  to  10  wave- 

lengths long  (reference  1).  Success  in  this  application  motivated  a 
study  of  multiple  wire  configurations  and  a test  of  application  moti- 
vated to  wire  grid  problems  (reference  2).  Initially,  studies  were 
^ restricted  to  problems  involving  100  unknowns  or  less  for  economic  reasons; 

exact  solutions  computed  by  standard  methods  were  compared  to  the  itera- 
tive solutions.  The  relative  efficiency  of  the  BMI  technique  was  suffi- 
ciently high  to  justify  construction  of  a computer  code  for  performing 
, the  out-of-core  storage  manipulations  required  in  the  solution  of  large 

problems.  During  construction  of  this  code,  the  numerical  properties  of 
the  technique  were  studied  and  convergence  measures  were  investigated 
(reference  3).  The  combination  of  the  use  of  symmetries  and  the  BMI 
study,  wire  grid  models  of  conducting  bodies  of  revolution  were  investi- 
gated. Model  parameters  were  varied,  and  computed  results  were  compared 
to  exact  theoretical  results.  Consistently  good  agreement  was  obtained, 
and  the  resulting  model  criteria  were  used  in  calculations  for  models  of 
objects  of  varied  shape  and  size.  Problems  with  up  to  1000  unknowns 
were  studied  (reference  5).  One  summary  paper  was  published  following 
the  study  of  single  wires  (reference  6),  another  was  published  following 
the  study  of  multiple  wires  (reference  7)  and  a third  was  published 
after  large  problems  were  studied  (reference  3). 

In  parallel  with  these  studies,  a GEMACS  computer  code  was 
undergoing  initial  development.  The  long  term  goal  is  to  incorporate  a 
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wide  range  of  methods  for  solving  electromagnetics  problems  in  a single 
user-oriented  package. 

2.  Theory 

The  method  of  moments  is  a formalism  for  reducing  an  integral 
equation  to  a set  of  linear  simultaneous  equations 


AX  = b, 


(Eq.  64) 


where  A is  the  complex  NxN  impedance  matrix,  X is  the  column  vector  of 
complex  coefficients  in  the  current  expansion,  and  b is  the  complex 
excitation  column  vector.  A variety  of  choices  for  the  integral  equa- 
tion, expansion  functions,  and  weighting  functions  are  in  use.  It  is 
assumed  that  the  combination  chosen  leads  to  an  unsymmetric  matrix  A. 

If  N is  sufficiently  small,  equation  64  can  be  solved  without 
using  peripheral  device  storage.  Comparative  efficiencies  of  different 
solution  algorithms  can  then  be  predicted  from  the  required  number  of 
complex  mo's  (multiplicative  operations).  If  N is  large,  other  factors 
must  be  considered  in  determining  relative  efficiencies.  These  are 
discussed  in  subsection  7. 

The  most  efficient  general  method  for  solving  linear  simul- 
taneous equations  is  to  decompose  the  matrix  into  a product  of  lower  and 
upper  triangular  matrices  using  Gaussian  elimination  (reference  9).  This 
requires  approximately  N^/3  mo's.  Solution  in  fewer  operations  requires 
some  special  feature  of  the  equations 

For  thin-wire  moments  problems  of  sufficient  size,  such  a 
special  feature  is  available.  The  matrix  elements  correspond  to  inter- 
actions between  wire  segments.  The  interactions  decrease  with  increasing 
distance  between  the  segments.  A detailed  in  subsection  C.5,  the  segments 
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can  be  numbered  such  that  the  difference  between  segment  numbers  for  all 
close-neighboring  segment  pairs  is  small  compared  to  N.  The  largest 
matrix  elements  can  then  be  kept  close  to  the  principal  diagonal  of  A. 
The  matrix  is  separated  into 

A = L + B + U,  (Eq.  65) 

where  B is  a banded  matrix  with  upper  and  lower  bandwidths  M (numbers  of 
minor  diagonals),  L is  the  triangular  portion  of  A below  B,  and  U is  the 
triangular  portion  of  A above  B.  Equation  64  can  be  written  as 


BX  = b - (L+U) X 


(Eq.  66) 


An  iterative  scheme  is  then 


BX. 


j+l  - 0 - 


(Eq.  67) 


where  X^  denotes  an  approximation  to  the  solution  at  the  iteration, 

and  X.  , denotes  an  approximation  to  the  solution  at  the  next  iteration, 
j + l 

Some  starting  value  Xq  is  chosen,  and  X^  is  computed  from  equation  67- 
Then  Xj  is  entered  on  the  right-hand  side,  and  X^  is  computed.  If  the 
sequence  converges,  an  approximate  solution  of  equation  64  is  obtained. 

Equation  67  must  be  solved  at  each  iteration.  The  cost  is 
minimized  by  decomposing  B into  a product  of  lower  and  upper  triangular 
banded  matrices,  B^  and  B^.  Equation  67  is  then  solved  by  forward 
elimination  in 


B,  Z.  = b - (L+U)  X . 
L J J 


(Eq.  68) 


followed  by  back  substitution  in 


Vj+i  a zj 


(Eq.  69) 
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Decomposition  of  B i s similar  to  full  matrix  decomposition,  except  that 

the  cost  can  be  much  less.  The  cost  depends  on  the  pivoting  strategy 

used.  Full  pivoting  destroys  bandedness,  and  and  B^  would  not  be 

banded.  Partial  pivoting  doubles  the  bandwidth  of  either  B^  or  B^. 

Pivoting  on  the  principal  diagonal  elements  of  B retains  the  bandwidth 

of  M for  both  B.  and  B ..  The  latter  is  the  least  expensive  in  terms  of 
l u 2 3 

mo's,  requiring  about  NM  - 2M  /3.  Pivoting  on  the  diagonal  elements 

increases  the  risk  of  large  errors  due  to  accumulating  rounding  errors 

during  decomposition.  This  subject  is  discussed  in  subsection  C.4. 

Assuming  that  B is  decomposed  by  pivoting  on  the  principal 

diagonal  elements,  equation  67  can  be  solved  at  a cost  of  about  N mo's 

for  each  iteration.  Assuming  that  K iterations  are  required  for  conver- 

2 3 

gence,  the  total  cost  for  the  iterative  solution  process  is  NM  - 2M  / 3 
2 

+ KN  mo's.  Based  on  the  number  of  mo's  required,  the  efficiency  g of 
the  banded  matrix  iterative  method  relative  to  the  best  general  method 


g = N3  [3  (NM2  - 2M3/3  + KN2)]’’ 


(Eq.  70) 


This  quantity  is  discussed  in  subsection  C.6. 

Theoretically,  convergence  of  the  sequence  is  assured  if  the 
spectral  radius  (magnitude  of  the  largest  eigenvalue)  of  B (L+U)  is 
less  than  one.  If  the  spectral  radius  is  greater  than  one,  the  sequence 
must  eventually  diverge.  If  it  is  only  slightly  greater  than  one,  the 
sequence  may  initially  converge  and  then  diverge.  This  behavior  is 
called  pseudoconvergence.  It  has  been  observed  in  this  research  (reference 
5).  The  bf’st  approximate  solution  obtained  during  pseudoconvergence  may 
be  sufficiently  accurate  for  some  purposes,  depending  on  the  quantity  or 
interest  and  the  percent  error. 
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For  the  iterative  method  to  be  efficient,  a reasonable  con- 
vergence rate  is  necessary.  If  the  choice  of  M and  the  segment  numbering 
scheme  results  in  such  a rate,  it  might  be  expected  that  the  solution  of 

| 

\ BX 1 = b (Eq.  71) 

t , 

would  yield  a reasonable  approximation  to  the  exact  solution  Xg.  This 

is  verified  in  practice.  As  a consequence,  the  zero  vector  can  be  used 

for  X , and  the  multiplication  (L+ll)X  need  not  be  performed.  The  cost 
° 2 ° 

of  obtaining  Xj  is  then  less  than  N (assuming  B has  been  decomposed). 

If  M<<N,  most  of  the  cost  of  one  iteration  is  saved.  This  is  signifi- 
cant only  for  small  problems.  The  important  point  is  that  if  the  con- 
vergence rate  is  high  enough  for  efficient  solution,  X^  is  a reasonable 
approximation  to  Xg  if  the  zero  vector  is  used  for  Xq.  Hence,  no  physical 
arguments  or  expensive  calculations  for  obtaining  a starting  value  are 
requi red. 

3 . Convergence  Measures 

The  iterative  solution  process  is  terminated  by  a numerical  CC 
(convergence  criterion)  on  some  convergence  measure.  In  the  early 
st'inies  of  small  problems  (references  1 and  2)  the  measure  used  was 
the  RE  (relative  error), 


REj  = [ (Xj  - Xe)+(Xj  - Xe)/XelXe]1/2,  (Eq.  72) 

where  (t)  denotes  the  complex  conjugate  transpose.  The  exact  solution 
Xg  was  obtained  by  full  matrix  decomposition  using  Gaussian  elimination. 
The  computer  CP  (central  processor)  times  for  example  problems  were 
recorded  both  for  the  exact  solution  process  and  for  the  BM I process. 
Efficiencies  computed  from  CP  times  and  from  numbers  of  mo's  (equation 
70)  were  similar,  verifying  that  the  number  of  mo's  is  an  adequate  mea- 
sure of  efficiency  for  solution  processes  if  no  out-of-core  man i pu 1 a t ions 
are  required. 


Due  to  the  cost  of  computing  for  large  problems,  the  RE  was 
available  as  a convergence  measure  only  during  the  research  phase  on 
small  problems.  An  alternative  measure  called  the  BCRE  (boundary  con- 
dition relative  error) 


BCREj  = [(AX.  - b) T (AX j 


b)/b+b]'/2 


(Eq.  73) 


was  proposed  (reference  1).  A study  of  the  relationship  between  the  BCRE 
and  the  RE  for  a variety  of  small  problems  showed  that  the  BCRE  is  not 
always  a reliable  measure  of  convergence  (reference  2).  This  lack  of  reli- 
ability can  be  traced  to  the  stability  of  the  equations  for  a given  prob- 
lem (reference  3)-  This  subject  is  discussed  in  subsection  C.b. 

An  alternative  measure  of  convergence  is  available  (reference  3) 
The  IRE  (iterative  relative  error)  is  defined  by 


IREj  = [(Xj  - Xj_|)+(Xj 


X.  , )/x.+x.] ,/2. 
J-)  J J 


(Eq.  7b) 


The  IRE  is  a measure  of  the  relative  change  in  the  sequence  of  approx- 
imate solutions  from  one  iteration  to  the  next.  It  is  an  adequate 
measure  of  convergence. 

The  sequence  of  values  of  the  RE  or  the  IRE  has  been  found  to 
be  approximated  by  a simple  exponential  function  (reference  3).  As  a 
result,  the  IRE  can  be  approximately  predicted  at  any  iteration  from 
the  values  at  the  previous  two  iterations; 


IRE.  = Pe_(i‘J 

1 P = 1 RE j (IRE 

IRE^  . = Pe_(i’  (J-1  ^ 
J-' 

0 

I UEJ  + | - Pe- 

Q = in ( 1 RE j _ 1 / 1 RE j ) 

© 

l*EJ  + | = ORE 

:j-r 


(Eq.  75) 
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The  smallness  criterion  for  convergence  can  then  be  placed  either  on  the 
IRE  or  on  the  predicted  IRE  at  the  next  iteration. 

Statistics  gathered  for  small  problems  showed  that  the  pre- 
dicted value  of  the  IRE  is  also  a good  prediction  of  the  RE,  at  least 

* for  cases  of  rapid  convergence  (reference  3).  This  leads  to  a definition 
of  the  ?RE  (predicted  relative  error) 

i 

• 2 

PREj  = ( I REj ) / IREj_ j • (Eg.  76) 

The  PRE  (Equation  75,  number  6)  is  an  adequate  measure  of  convergence. 

The  flexibility  of  the  GEMACS  code  allows  the  user  to  specify  his 
own  convergence  measure.  Any  quantity  that  is  readily  (inexpensively) 
computed  can  be  used.  Since  some  quantities  such  as  the  BCRE  may  not  be 
adequate  due  to  instability  of  the  equations,  new  measures  should  be 

' used  wi th  caution. 

4.  Stab! 1 i ty 

The  subjects  in  this  section  have  been  investigated  in  great 
detail  in  the  last  two  decades.  The  purpose  of  this  section  is  to 
introduce,  in  a nonrigorous  manner,  those  definitions  and  concepts 
which  were  useful  in  this  research.  The  material  and  notation  closely 
follow  reference  9-  Most  quantities  in  this  report  are  complex,  whereas 
those  in  reference  9 are  real. 

The  euclidean  length  or  norm  of  a vector  X in  complex  Ni- 
di mens  ional  space  is  defined  as 

I | X | | = (X+X)'/2  (Eq.  77) 

Other  vector  norms  exist,  but  will  not  be  used  here.  The  norm  of  a 
complex  matrix  A with  N rows  and  N columns  is  defined  as 
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INI  = max  J lyT  tf,/  l (Eq.  78) 

X*0  |X  I V» 

'X 

where  0 denotes  the  zero  vector.  An  alternative  definition  is 


I I AX | | , | | X | |=1 


(Eq.  79) 


This  can  be  interpreted  to  mean  that  if  the  unit  sphere  in  the  space  of 
X is  mapped  by  AX  into  the  space  of  b,  ||a||  is  the  longest  vector  that 
will  be  obta i ned . 

Similarly,  the  norm  of  the  inverse  of  A can  be  defined  by 


| | A '||  = max 


(Eq.  80) 


or  by 


- 1 I A~ 1 b | | , ||b||=l 


(Eq.  81) 


Thus,  if  the  unit  sphere  in  the  space  of  b is  mapped  into  the  space  of  X 
by  A 'b,  | |A  ' | | is  the  longest  vector  that  will  be  obtained. 

The  condition  number  of  A,  denoted  by  cond (A) , is  defined  as 


cond(A)  = 


(Eq.  82) 


This  number  is  not  easily  computed.  It  is  conceptually  valuable  in  a 
variety  of  ways.  If  the  vector  b is  subject  to  the  uncertainty  <5b,  the 
uncertainty  6X  in  X is  bounded  by 


6X  l/nN  6b 

Tx  - cond(A)  • -TV 


(Eq.  83) 


The  relative  error  in  the  solution  is  then  bounded  by  the  product  of  the 
relative  error  in  b and  cond(A).  If  the  matrix  elements  in  A are  subject 
to  uncertainty  6A,  a similar  bound  is  obtained.  In  each  case,  a large 
value  of  cond  (A)  is  a warning  that  the  solution  may  be  highly  sensitive 
to  small  changes  in  A or  b. 
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Another  useful  bound  involves  the  ratio  of  the  RE  to  the 
relative  residual,  which  is  the  BCRE  used  in  this  work: 

1 RF 

i Bcs^icond<A,•  w 

Since  cond(A)  >_  1,  this  ratio  can  be  very  small  or  very  large.  For  a 
variety  of  thin-wire  moments  problems,  this  ratio  was  found  (reference  2) 
to  be  between  1 and  10,  with  an  average  of  about  3.  It  is  obvious  that  a 
much  wider  range  of  values  is  possible.  If  cond(A)  is  large,  A is  said 
to  be  ill-conditioned.  One  consequence  of  ill-conditioned  matrices  is 
that  the  BCRE  may  be  quite  small  when  the  RE  is  large.  The  BCRE  may 
then  be  an  adequate  measure  of  convergence  in  the  banded  matrix  iter- 
ation method  only  if  cond (A)  is  near  unity. 

A measure  of  stability  that  has  been  used  in  the  method  of  moments 
(reference  10)  and  that  has  properties  similar  to  those  of  the  condition 
number,  is  the  pivot  ratio.  This  quantity  is  easily  calculated  during 
decomposition  by  Gaussian  elimination.  It  can  be  defined  as  the  ratio 
of  the  magnitudes  of  the  first  and  last  pivot  elements  (reference  10),  or 
as  the  ratio  of  the  largest  to  the  smallest  of  the  magnitudes  of  the  pivot 
elements.  Only  the  order  of  magnitude  of  the  pivot  ratio  is  significant. 

The  choice  of  definitions  is  not  important  when  compared  to  the  effect 
of  the  choice  of  pivoting  strategies.  While  pivoting  on  the  largest 
element  of  a row  or  column  will  reduce  the  cumulative  effect  of  rounding 
errors,  it  also  reduces  the  pivot  ratio.  Hence,  pivoting  on  the  diagonal 
elements  should  produce  a larger  pivot  ratio  for  ill-conditioned  matrices. 

In  the  banded  matrix  iterative  scheme,  the  banded  matrix  is 
decomposed  us!ng  Gaussian  elimination.  The  decomposition  is  accomplished 
by  pivoting  on  the  diagonal  elements.  (Pivoting  on  other  elements  would 
increase  the  bandwidth  and  result  in  loss  of  efficiency.)  The  pivot 
ratio  for  the  banded  matrix  is  then  available  as  a measure  of  ill-con- 
ditioning of  B.  Unfortunately,  no  method  has  been  found  to  relate 
either  the  pivot  ratio  to  cond(B)  or  cond(A)  to  cond(B).  It  is  possible 
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that  either  A or  B could  be  ill-conditioned  while  the  other  is  not. 

This  would  probably  depend  on  the  physical  problem,  the  choice  of  segment 
numbering  schemes,  and  the  choice  of  bandwidths.  Certainly  cond(B) 
should  approach  cond (A)  as  M approaches  N. 

5.  Segment  Numbering 

Two  types  of  segment  numbering  problems  occur.  For  the  iter- 
ative scheme  to  be  efficient,  the  segments  must  be  numbered  so  as  to 
keep  the  large  matrix  elements  in  the  band.  The  best  choice  of  num- 
bering for  this  purpose  is  obvious  in  some  cases  and  not  in  others. 

More  than  one  choice  may  be  apparent.  The  same  logical  process  can  be 
used  to  number  segments  for  a variety  of  problems,  as  will  be  shown. 

The  other  problem  in  numbering  is  related  to  the  requirements  for  input 
data  for  the  geometry  processor  in  any  moments  code.  Generally  speaking, 
the  code  automat i ca 1 1 y numbers  segments  in  sequence  along  each  wire,  in 
the  order  in  which  wire  data  are  entered.  For  the  iterative  method, 
this  numbering  is  adequate  in  some  cases  and  not  in  others.  If  some 
other  numbering  is  necessary,  the  obvious  method  is  to  enter  the  wire 
data  in  the  order  desired.  For  most  codes,  this  could  require  a tedious 
process  of  entering  each  segment  as  a single  wire  and  supplying  connection 
data.  An  easier  way  will  be  shown  that  involves  segment  renumbering. 

a . Numbering  for  the  Iterative  Method 

The  logic  for  the  numbering  scheme  that  is  adequate  for 
many  problems  is  most  easily  explained  for  flat  objects  (all  wires  in  a 
plane).  The  basic  idea  is  to  superimpose  a set  of  narrow  parallel 
strips  on  the  object  of  interest.  Figure  3 shows  the  strips  as  separated 
by  dashed  lines,  with  a lumpy  appearing  object  in  outline.  Assuming 
that  the  object  is  a wire  loop,  numbering  is  as  shown  in  the  figure. 
Starting  at  one  extremity,  numbering  proceeds  from  left  to  right  until 
all  segments  in  the  first  strip  are  numbered.  Moving  to  the  next  strip, 
numbering  again  proceeds  from  left  to  right.  This  continues  until  all 
segments  are  numbered.  In  the  figure,  it  is  assumed  that  the  segment 
end  occurs  at  intersections  with  the  dashed  lines.  This  need  not  be  the 
case.  Strip  widths  can  also  vary.  This  depends  on  the  geometry  of  the 
object . 


^5 


r 


Section  E of  reference  2 includes  a discussion  of  this 
numbering  scheme  as  it  applies  to  example  problems.  Generally,  any  flat 
' object  should  be  oriented  so  that  the  strips  run  across  the  narrow 

dimension  of  the  object.  This  keeps  the  number  difference  between  seg- 
ments in  adjacent  strips  as  small  as  possible. 

The  strip  numbering  scheme  has  a direct  application  for 
‘ some  three-dimensional  objects.  The  simplest  is  a cylindrical  grid. 

Numbering  a rectangular  grid  by  the  strip  method  and  rolling  the  grid 
into  a cylinder  (about  the  correct  axis)  results  in  a helix-like  number- 
ing scheme.  This  is  appropriate  when  the  cylinder  length  is  greater 
than  or  roughly  equal  to  its  circumference.  If  the  cylinder  has  end 
caps,  numbering  should  proceed  from  one  end  cap  center  to  the  other. 

The  pattern  is  then  spiral,  helical,  and  spiral.  However,  if  the 
cylinder  is  short  compared  to  its  ci rcumference,  the  spiral  scheme  may 
t not  give  an  efficiency  as  high  as  that  obtained  by  orienting  the  cylinder 

axis  normal  to  the  plane  of  the  strips. 

A sphere  or  cube  can  be  handled  much  the  same  as  the 
cylinder  with  end  caps.  For  elongated  shapes,  numbering  should  start  at 
• one  extremi ty . 

A similar  method  for  segment  numbering  that  can  be  auto- 
mated is  discussed  in  reference  A.  It  is  called  geometric  cell  divi- 
sion. Figure  A shows  a set  of  parallel  strips  normal  to  a direction 
vector  d.  For  planar  objects,  is  was  noted  that  the  set  of  strips  could 
be  superimposed  on  the  wire  object;  segment  numbering  would  proceed 
along  each  strip  starting  at  one  extremity  of  the  object.  For  an  irregu- 
lar object,  the  choice  of  direction  for  d is  not  obvious. 

Figure  5 shows  an  irregular  planar  object  with  four 
possible  directions  for  d.  If  the  small  appendages  were  absent,  d^  would 
be  the  obvious  choice.  Numbering  would  then  be  in  strips  across  the 
narrow  dimension  of  the  large  rectangle.  With  the  appendages  present, 
directions  d^  and  d^  are  not  advisable.  In  some  strips,  numbering  would 
proceed  across  the  large  rectangle  and  along  the  length  of  one  projection 
before  advancing  to  the  next  strip.  This  can  be  alleviated  somewhat 
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by  using  direction  d2>  so  that  no  strip  runs  through  the  length  of 
the  longer  projection.  The  other  choice  is  d^,  and  it  is  clearly  the 
best.  No  strip  runs  both  across  the  large  rectangle  and  along  a projec- 
tion, and  the  strips  that  cut  through  both  projections  include  few  wire 
segments. 

*d 


Figure  A.  A Set  of  Parallel  Strips  Normal  to  a Direction  Vector  d 

After  studying  several  examples  of  this  sort,  a consistent 
result  emerges.  The  best  choice  of  orientation  apparently  lies  along 
(or  at  some  small  angle  to)  a principal  inertial  axis  of  equal  point 
masses  located  at  the  wire  segment  centers.  The  best  principal  axis  is 
always  the  one  about  which  the  moment  of  inertia  is  least.  This  scheme 
will  be  called  the  PASS  (Principal  Axis  Slicing  System).  No  proof  has 
been  found  that  this  orientation  for  slicing  is  the  best  choice,  and 
its  practicality  would  have  to  be  exhibited  in  practice.  The  method  can 
be  extended  to  three  dimensions  as  follows.  Let  the  principal  axis 
coordinates  be  denoted  X,  Y,  Z,  with  the  least  inertia  about  the  X-axis 
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and  the  greatest  about  the  Z-axis.  Divide  the  object  by  a set  of  planes 
normal  to  the  X-axis,  with  the  separation  of  the  planes  somewhat  less 
than  the  wire  segment  lengths.  The  space  between  adjacent  planes  is 
divided  by  planes  normal  to  the  Y-axis,  leading  to  a set  of  parallel 
tubes  of  square  cross  section  in  each  planar  section  normal  to  the  X- 
axis.  Numbering  proceeds  along  one  tube  (parallel  to  the  Z-axis)  starting 
at  one  edge  of  the  section  of  the  object.  It  proceeds  from  tube  to  tube 
l across  the  object,  and  then  to  the  next  planar  section.  The  first 

planar  section  is  always  at  one  extremity  of  the  object.  Figure  6 
shows  the  numbering  obtained  by  this  method  for  a rectangular  object 
modeled  with  a regular  wire  grid. 

For  a problem  such  as  a large  cylinder  with  end  caps  or 
other  such  flat-ended  object,  the  method  of  geometric  cell  division  is 
probably  inferior  to  the  spiral-helical-spiral  numbering  method.  The 
former  results  in  numbering  across  the  end  in  strips  before  proceeding 
, along  the  object,  so  that  a rather  large  difference  in  segment  numbers 

occurs  near  the  edge  where  numbering  begins.  The  limited  practical 
experience  with  large  fat  objects  precludes  any  definite  statement  con- 
cerning the  best  approach. 

• b.  Segment  Renumbering 

As  mentioned  earlier,  the  geometry  processing  portions  of 
most  programs  provide  segment  numbering  that  is  sequential  along  each 
wire  in  the  order  in  which  wire  data  are  entered.  This  segment  number- 
ing is  adequate  for  the  banded  matrix  method  for  some  problems.  If  it 
is  not,  it  can  still  be  used  initially  to  simplify  the  model  input  data. 
The  segments  can  then  be  renumbered  in  any  sequence  desired.  Details 
and  examples  of  the  renumbering  process  as  it  was  used  in  the  research 
phase  are  contained  in  Section  E of  reference  2.  The  GEMACS  code  permits 
the  user  to  number  the  segments  in  any  order  desired. 
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Figure  6.  Part  of  the  Segment  Numbering  Obtained  with  the  Principal 
Axis  Slicing  System  for  a Rectangular  Pa ra 1 1 e 1 p i ped 
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6.  Choice  of  Bandwidth 

The  choice  of  bandwidth  M depends  both  on  the  problem  and  on 
the  numbering  scheme  used.  The  ideal  choice  yields  a minimum  solution 
cost.  An  understanding  of  the  formula  for  efficiency  in  equation  70  is 
helpful  in  making  this  choice. 


Efficiency  Characteristics 
The  efficiency  can  be  written 


] , (Eq.  85) 

f2(3-2f)  + 3K/N 


where  f is  M/N.  A useful  characteristic  of  this  efficiency  is  obtained 
by  setting  K equal  to  zero.  The  upper  limit  on  g is  then 


This  function  is  shown  in  figure  7 as  a solid  line,  treating  f as  a 
continuous  variable.  The  actual  efficiency  always  falls  on  or  below 
this  line. 

For  any  given  problem,  with  a particular  choice  of  segment 
numbering  and  convergence  criterion,  the  actual  efficiency  as  a function 
of  the  bandwidth  is  as  follows.  At  M equal  to  N,  the  solution  will 
always  satisfy  the  CC  (convergence  criterion).  (When  the  relative  error  is 
not  available  as  a convergence  measure,  at  least  one  iteration  is  necessary 
to  test  for  convergence.  If  N is  large,  one  iteration  is  of  little 
cost,  and  will  be  ignored  in  this  discussion.)  As  M is  decreased,  the 
error  in  Xj  gradually  increases.  Experience  shows  that  this  error  may 
not  increase  monoton i ca 1 1 y , but  may  have  minor  oscillations  superimposed 
(reference  1).  Eventually,  a bandwidth  in  the  breaking  region  is  reached. 
This  region  is  a narrow  range  in  bandwidth  where  the  error  in  Xj  is  just 
above  or  just  below  the  CC.  A breaking  point  can  be  defined  as  the  largest 
bandwidth  at  which  the  error  in  X^  is  greater  than  the  CC.  At  this 
bandwidth,  at  least  one  iteration  is  required  to  reduce  the  error  to 
the  CC.  In  most  cases,  the  actual  efficiency  will  then  depart  from 
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f = 'M/N 

Figure  7.  The  Curve  of  g. ...  Versus  f and  a Fictitious 
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Example  of  Actual  Efficiencies 
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the  gLjM  curve.  Because  the  savings  in  cost  of  decomposing  B at  a 
smaller  bandwidth  is  greater  than  the  cost  of  a few  iterations  in  the 
breaking  region,  the  efficiency  generally  continues  to  rise  as  M decreases. 
The  rate  of  increase  is  smaller  than  that  of  g. lu,  however. 

L I n 

^ In  some  cases  the  actual  efficiency  can  depart  from  §L|M 

:>  at  the  breaking  point,  then  return  briefly  due  to  one  of  the  minor 

oscillations  mentioned  earlier,  and  then  depart  again.  This  behavior  was 
shown  by  a fictitious  example  in  figure  7. 

For  any  particular  CC,  the  actual  efficiency  will  rise 
somewhat  from  the  value  at  the  breaking  point  as  M decreases,  then  reach 
a peak  and  begin  to  decline.  The  discrete  nature  of  the  process  neces- 

% 

sarily  means  that  minor  up  and  down  behavior  over  narrow  ranges  of  band- 
widths  will  occur.  This  behavior  is  highly  dependent  on  the  problem  and 
the  segment  numbering  scheme. 

Figure  8 shows  the  curve  of  g.  ...  and  a set  of  dashed 
curves  obtained  from  equation  70,  using  the  data  from  example  problem  9 
in  reference  2.  (The  convergence  measure  was  the  RE.  Points  along  the 
dashed  curves  were  computed  from  the  data  in  Table  A5  of  reference  2, 
reproduced  here  as  table  1.  Points  for  the  curve  at  a CC  of  0.1  percent 
were  approximated  using  the  extended  exponential  function.) 

Consider  the  curve  for  a CC  of  10  percent.  At  f equal  to 
17/80  or  0.21,  the  RE  for  was  6.13  percent.  This  satisfies  the 
CC  of  10  percent,  so  no  iterations  were  required.  Then  either  equation 
70  or  equation  86  yields  an  efficiency  of  8.7-  This  is  the  highest 
efficiency  that  can  be  obtained  for  this  problem  at  a bandwidth  of  17; 
relaxing  the  CC  will  not  affect  the  efficiency.  Tightening  it  to  some- 
thing less  than  6.13  percent  will  result  in  at  least  one  iteration  being 
performed,  with  a loss  in  efficiency. 

The  location  of  the  peak  of  the  efficiency  curve  clearly 
shifts  toward  smaller  bandwidths  as  the  CC  is  relaxed.  The  width  of  the 
curve  near  the  peak  narrows  as  the  CC  is  relaxed.  The  rate  of  decrease 
in  efficiency  is  less  to  the  right  of  the  peaks. 
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TABLE  1.  RELATIVE  ERROR  (%)  FOR  EXAMPLE  9,  REFERENCE  2 (VERTICAL 


HALF-RHOMBIC  ANTENNA,  N=80) 


ITERATIONS 

4 

5 

BANDWI DTH 

6 7 

12 

17 

1 

67.67 

53-91 

40.47 

38.91 

6.13 

2 

46.09 

28.89 

16.  14 

15- 12 

m 

.59 

3 

31.62 

15-64 

6.53 

5-94 

.18 

A 

21.69 

8.46 

2.64 

2.33 

5 

14.87 

4.57 

1.06 

• 91 

6 

10.20 

2.47 

• 43 

7 

7.00 

1.34 

8 

4.80 

• 72 

9 

3-29 

10 

2.26 

1 1 

1.55 

12 

1 .06 

13 

• 73 
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If  the  program  user  imposes  a tight  CC  (about  1 percent 
or  smaller  for  small  problems),  he  can  afford  to  pick  M such  that  he 
misses  the  peak  efficiency  by  quite  a bit  to  either  side  because  the 
peak  is  wide.  At  larger  values  of  the  CC,  he  will  tend  to  decrease  M to 
try  to  stay  near  the  peak,  but  he  may  be  increasingly  conservative  to 
avoid  a point  too  much  to  the  left  of  the  peak.  This  is  good  in  a 
sense,  because  use  of  larger  bandwidths  provides  lower  initial  errors 
and  higher  convergence  rates,  resulting  in  better  performance  of  the  PRE 
convergence  measure  (see  subsection  C.3).  If  the  user  picks  a bandwidth 
too  much  to  the  left  of  the  peak,  slow  convergence  and  low  efficiency 
will  occur.  Restarting  with  a larger  bandwidth  might  provide  a higher 
overall  efficiency  even  with  the  cost  of  an  aborted  run  included, 
b . Bandwidth  Estimates  for  Long,  Thin  Objects 

The  shapes  and  locations  of  peaks  of  the  efficiency 

curves  vary  considerably  from  problem  to  problem.  No  method  for  accurate 

prediction  of  the  bandwidth  for  peak  efficiency  has  been  found.  A 

survey  of  example  small  problems  did  yield  a trend  when  the  CC  is  1 

percent.  It  was  found  that  the  bandwidth  for  efficient  solution  could 

be  estimated  from  structure  dimensions  with  fair  reliability.  A bandwidth 

M corresponds  to  a distance  R..  within  which  all  interactions  are  to  be 

n 

kept  in  the  band.  (For  segments  of  length  0.1  X numbered  sequentially 
along  a straight  wire,  a bandwidth  of  M corresponds  to  a distance  of  H x 
0.1  X.)  The  distance  in  wavelengths  for  bandwidth  estimation  is  shown 
in  figure  9 as  a function  of  the  object  length  in  wavelenqths  for 
objects  having  a dominant  dimension  L.  The  vertical  bars  show  the 
uncertainty  in  value  for  several  problems.  The  linear  trend  is  obvious. 

This  trend  was  obtained  from  studies  of  problems  with  100 
unknowns  or  less.  They  were  primarily  problems  involving  long  wires  or 
wire  arrays.  Studies  of  larger  problems  did  not  include  bandwidth  as  a 
parameter.  The  initial  bandwidth  was  generally  selected  from  figure  9. 

For  planar  or  near-planar  objects  and  for  thin  cylinders,  rapid 
convergence  was  fairly  consistently  obtained.  Fat  cylinders,  spheres, 
and  objects  where  L is  less  than  X require  larger  bandwidths  (reference  5) • 
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Figure  9.  The  Linear  Trend  Between  the  Bandwidth 
Estimate  and  the  Object  Length  L 
(R  = Distance  Corresponding  to  Bandwidth) 
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The  parameter  L is  usually  the  largest  geometric  dimen- 
sion; however,  for  complex  objects,  the  choice  or  value  of  L may  depend 
on  the  excitation  type  and  orientation  (this  is  shown  by  an  example 
in  subsection  D.3).  Use  of  figure  9 for  these  objects  may  lead  to 
pseudo-convergence  for  cases  in  which  there  are  large  (or  small  coherent) 
interactions  not  included  in  the  band. 

7 . Computer  Timing  Requirements 

For  small  problems,  the  primary  computer  cost  is  the  CP  (central 
processor)  time  required  to  compute  matrix  elements.  Use  of  the  BMI 
solution  process  in  the  place  of  full  matrix  decomposition  by  Gaussian 
elimination  is  primarily  of  academic  interest  for  such  problems.  The  CP 
time  for  computing  matrix  elements  increases  as  the  square  of  N for  most 
programs.  The  solution  time  increases  as  the  cube  of  N for  full  matrix 
decomposition.  For  sufficiently  large  N,  the  solution  time  by  this 
method  will  eventually  dominate.  The  iterative  process  is  a method  for 
reducing  the  solution  cost.  The  total  cost  for  generating  and  solving 
large  problems  is  sufficient  to  justify  optimization  of  the  code  for  a 
given  computing  facility. 

For  problems  with  N in  the  range  below  1000,  three  primary 
factors  contribute  to  the  computer  costs;  the  CP  time  required  to  compute 
or  fill  the  impedance  matrix,  the  CP  time  required  to  perform  the  multi- 
plicative operations  in  solving  the  equations,  and  the  PP  (peripheral 
processor)  time  required  for  out-of-core  manipulations.  The  PP  time 
includes  input/output  time  and  considerable  bookkeeping  time.  For  large 
problems,  the  PP  time  is  almost  entirely  associated  with  the  solution 
process.  A cost  comparison  of  solution  methods  should  include  both  CP 
and  PP  times.  (The  best  comparison  is  actually  based  on  dollar  costs, 
but  no  two  installations  use  the  same  cost  algorithm  to  charge  for  total 
computer  resources  used.)  The  PP  time  and  solution  CP  time  are  extremely 
dependent  on  the  machine,  the  computing  system,  and  the  particular 
computer  code  being  used.  The  solution  CP  time  also  depends  on  the 
compiler  used.  Consequently,  the  efficiency  based  on  the  number  of 
multiplicative  operations  is  used  as  a simple  assessment  of  overall 
ef f i c i ency . 


The  CP  times  required  for  solution  on  the  CDC  6600  computer 
(FTN  k.2  compiler)  are  shown  in  figure  10  versus  the  number  of  multi- 
plicative operations  required  to  solve  example  problems  using  the  itera- 
tive method.  The  linear  dependence  over  such  a wide  range  indicates  that 
the  efficiency  based  on  number  of  multiplicative  operations  is  a valid 
comparison  of  CP  costs  between  methods. 

The  PP  time  shows  a similar  linear  dependence  on  the  number  of 
multiplicative  operations,  except  that  the  PP  times  were  higher  by  a 
factor  of  about  15.  This  proportionality  factor  is  extremely  dependent 
on  the  amount  of  fast  access  core  available  for  matrix  elements  during 
the  decomposition  of  B.  (For  a given  amount  of  storage,  the  PP  time  for 
decomposition  of  B should  be  much  less  than  that  for  decomposition  of  the 
full  matrix  using  Gaussian  elimination.)  The  amount  of  available  core  can 
be  increased  by  a variety  of  methods  including  program  segmentation  or 
overlaying.  These  methods  were  not  used  during  this  investigation  and 
only  1A  columns  of  the  matrix  could  be  kept  in  core  for  a problem  with 
1000  unknowns.  As  a consequence,  the  PP  costs  were  the  highest  dollar 
cost  of  the  study. 

The  CP  times  required  by  the  modified  AMP  code  incorporating 
the  BMI  solution  technique  for  computing  matrix  elements  is  shown  in 
figure  11  versus  the  number  of  unknowns.  The  upper  curve  is  for  wire 
grid  problems  and  the  lower  is  for  wire  problems  with  no  multiple  wire 
junctions.  The  difference  is  due  to  a repeated  search  of  junction  con- 
nection data,  and  is  not  apparent  for  small  problems.  Most  of  this  cost 
difference  is  eliminated  in  the  GEMACS  code  by  replacing  the  search  with 
a circular  linked  list.  The  CP  times  for  large  problems  can  be  further 
reduced  by  using  less  expensive  algorithms  (such  as  the  Hertzian  dipole) 
to  compute  matrix  elements  corresponding  to  interactions  between  segments 
separated  by  about  a wavelength  or  more.  This  option  was  not  used  for 
the  example  problems.  Further  cost  reductions  may  be  obtained  by 
relaxing  the  accuracy  restrictions  in  computing  large  matrix  elements. 

These  elements  should  be  accurate  only  enough  to  match  the  overall  accuracy 
required  of  the  modeling  and  solution  process. 
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D.  SUMMARY  OF  RESULTS 


Numerical  results  were  obtained  for  a variety  of  thin-wire  problems 
(references  1 through  5).  Three  different  computer  programs  were  used. 

The  general  results  and  a few  examples  are  listed  here.  All  calculations 
were  performed  on  a CDC  6600  computer.  Dimensions  are  in  meters. 

1 . Single  Straight  Wires 

The  first  results  (reference  1)  were  obtained  using  the  moments 
formalism  of  Harrington  (reference  11).  That  is,  the  potential  integral 
equation  was  solved  using  pulse  expansion  functions  and  point  matching 
(collocation)  at  wire  segment  midpoints.  All  examples  were  single 
straight  wires,  either  antennas  or  broadside  scatterers,  with  L (length) 
to  X (wavelength)  ratios  ranging  f "om  2.5  to  10.  Length  to  D (diameter) 
ratios  were  either  7^*2  or  1000.  The  con 'ergence  measure  was  the  RE, 
and  the  convergence  criterion  was  i percent.  There  were  10  segments  per 
wavelength  used,  with  segment  numbering  sequential  along  the  wire.  It 
was  found  that: 

(1)  The  iterative  solution  converges  monoton ica 1 1 y at  bandwidths 
from  3 to  N.  (With  10  segments  per  wavelength,  a bandwidth  of 
3 means  that  interactions  at  distances  over  0.3A  are  excluded 
from  the  band.)  An  example  of  the  convergence  behavior  is 
shown  in  figure  12.  The  wire  is  a centerfed  antenna  with  i/X 
= 2.5,  L/D  = 1000,  and  N = 25.  The  figure  shows  the  RE  at 
each  iteration  for  a given  bandwidth.  Points  are  connected  by 
a line  for  clarity;  the  results  are  necessarily  discrete. 

(2)  Convergence  rates  depend  on  the  wire  radius,  as  do  errors  in  the 
initial  solution,  Xj.  Thinner  wires  result  in  better  efficiency. 

(3)  Relative  errors  for  are  larger  for  resonant  antennas  than 
for  nonresonant  antennas  for  the  same  bandwidth.  The  differ- 
ence is  greater  for  smaller  wire  radii. 

(A)  The  relative  error  for  does  not  decrease  smoothly  and 

montonically  with  increasing  bandwidth.  Periodic  variations 
occur  at  changes  in  bandwidth  that  correspond  to  half-wavelength 
changes  in  distance  along  the  wires,  as  exemplified  in  figure 
13. 
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(5)  As  indicated  in  subsection  C.3,  convergence  behavior  for  a 
given  bandwidth  can  be  approximated  with  varying  accuracy  by  the 
formula  RE(%)  100  e ^ where  i is  the  number  of  iterations 

and  Q is  a function  that  is  dependent  on  the  bandwidth.  The 
dependence  is  linear  in  a first  approximation,  with  an  oscilla- 
tion superimposed.  However,  the  slope  of  the  linear  part  is 
dependent  on  the  excitation,  length,  and  wire  radius. 

(6)  For  a length  to  wavelength  ratio  of  10,  iterative  solutions 
with  a relative  error  of  1 percent  were  obtained  with  an 
efficiency  of  about  7.  For  a relative  error  of  10  percent,  an 
additional  reduction  in  CP  time  up  to  a factor  of  2 was  obtained. 

(7)  The  bandwidth  resulting  in  peak  efficiency  varies  with  wire 
length  in  a roughly  linear  fashion.  At  10  segments  per  wave- 
length, the  optimum  bandwidth  is  about  4 of  5 for  3*  wires  and 
in  the  range  10-15  for  10A  wires. 

The  reader  is  referred  to  reference  1 for  details  and  to  6 and  7 for 
summary  comments. 

2 . Multiple  Straight  Wires 

The  second  set  of  results  (reference  2)  was  obtained  by  modi- 
fying program  WAMP  (reference  12).  That  program  is  based  on  the  Pockling- 
ton  integral  equation,  with  pulse  plus  sine  plus  cosine  expansion  func- 
tions and  collocation.  With  N wire  segments,  this  choice  of  expansion 
function  requires  3N  current  coefficients.  Imposition  of  "extended 
continuity  conditions"  at  adjacent  segment  midpoints  reduces  the  number 
of  unknowns  to  the  il  values  of  current  at  the  segment  midpoints.  Point 
continuity  in  the  current  is  not  obtained  at  junctions.  An  extended 
continuity  method  is  also  used  at  multiple  junctions,  but  the  method 
is  considered  poor  for  wire  grid  problems. 

Using  this  program  for  generating  the  equations,  various 
combinations  of  thin-wire  geometries  were  used  to  investigate  the 
capabilities  and  limitations  of  the  banded  matrix  method.  These 
i ncl uded : 
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(1)  One  straight  wire. 

(2)  Two  parallel  centerfed  antennas  at  varying  separations  and 
rad  i i . 

(3)  Two  col  linear  centerfed  antennas  at  large  and  smail  separations. 

(4)  A linear  array  of  parallel  dipoles. 

(5)  Square  and  circular  arrays  of  parallel  dipoles. 

(6)  A two-dimensional  array  of  short  dipoles. 

(7)  A hel ix  antenna . 

(8)  A vertical  half-rhombic  antenna  over  sea  water. 

(9)  A wire-gridded  rectangular  strip. 

(10)  A square  loop. 

(11)  A pair  of  crossed  wires. 

In  each  case,  N was  restricted  to  100  or  less  for  economic  reasons. 

For  many  of  the  examples,  the  sequence  of  solutions  converged 
uniformly  to  the  exact  solution.  Divergence  was  forced  for  some  examples 
by  using  a combination  of  element  numbering,  geometry,  and  bandwidth  so 
that  some  large  matrix  elements  were  not  contained  in  the  B band.  The 
general  approach  to  segment  numbering  for  the  BM I solution  technique  was 
developed  during  this  phase. 

Because  of  the  wide  variety  of  geometries  involved  in  these 
examples,  the  results  are  not  in  a form  to  be  easily  summarized.  The 
convergence  measure  was  the  RE.  At  a convergence  criterion  of  1 percent, 
efficiences  were  generally  around  5 to  10  and  up  to  23.  The  iterative 
process  was  interrupted  for  some  problems  when  the  RE  was  first  reduced 
to  less  than  30  percent.  The  far  field  pattern  was  computed  from  the 
approximate  solution.  The  iteration  was  then  continued  to  the  1 
percent  level,  and  the  far  field  pattern  was  again  computed.  Comparison 
of  the  far  fields  obtained  from  the  two  current  distributions  indicates 
that  a 30  percent  CC  may  be  adequate  for  far-field  pattern  determination. 

Example  12  of  reference  2 is  reproduced  here  to  exhibit  one 
of  the  problems  with  automating  the  segment  numbering  scheme  called  PASS 
(see  section  C.5).  Table  2 shows  the  description  of  the  problem;  a 


TABLE  2.  DESCRIPTION  OF  EXAMPLE  12  FROM  REFERENCE  2 


TYPE:  Square  loop 


DIMENSIONS:  2.5  meters  on  each  side 

Rw=  0.0015,  A = 1 


SEGMENTATION:  10  segments  per  meter,  total  of  100 


SEGMENT  NUMBERING: 


135  ...  49 


Excitation:  Segment  25  (center  of  bottom  edge) 

EXACT  SOLUTION  CP  TIME  (sec):  4.178 
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square  loop  antenna.  Dimensions  are  in  meters.  The  strip  numbering 
scheme  produces  the  segment  numbering  shown  if  the  strips  are  oriented 
properly.  (The  vector  d in  figure  A would  be  tilted  A5  degrees  to 
the  right.)  Table  3 shows  the  convergence  data  obtained.  Efficiencies 
are  shown  based  both  on  comparative  computer  central  processor  times  and 
on  comparative  numbers  of  multiplicative  operations. 

The  problem  with  using  an  automated  PASS  system  for  this 
problem  is  that  the  geometry  is  redunant.  The  principal  axes  can  lie 
either  parallel  to  the  sides  of  the  square  or  pass  through  its  corners. 

The  former  would  yield  a poor  numbering  scheme,  because  numbering  would 
be  sequential  along  one  edge,  then  alternate  between  opposite  sides.  An 
efficiency  greater  than  three  is  unlikely  with  such  numberinq. 

It  would  be  possible,  however,  to  automate  a slicing  system 
for  numbering  that  would  allow  the  user  to  specify  the  orientation  for 
the  vector  d. 

For  the  numbering  actually  used,  the  bandwidth  for  efficient 

solution  could  have  been  estimated  from  figure  9 using  an  object 

length  (corner  to  corner  circumferential  distance)  of  5A.  A distance  R., 

M 

of  about  0.7^  is  suggested  by  figure  9.  The  difference  in  segment 
numbers  for  segments  separated  by  0.7A  is  1**.  A bandwidth  of  lA  provided 
peak  efficiency. 

3 . Mod i f i ed  AMP 

The  results  documented  in  references  1 and  2 were  sufficiently 
encouraging  to  motivate  development  of  a computer  program  for  solving 
large  problems.  The  AMP  (Antenna  Model irg  Program)  was  selected  for 
modification  (reference  13).  The  AMP  cooe  uses  the  same  moments  formalism 
as  the  WAMP  code  except  for  continuity  schemes  at  junctions.  Reference  A 
discusses  the  modifications  and  the  resulting  code. 

A number  of  problems  of  intermediate  size  (N  <_  300)  was  investi- 
gated (reference  A).  The  most  important  of  these  were  the  wire  grid  models 
of  conducting  objects  of  various  shapes.  A parameter  study  (wire  radius, 
segment  length  and  associated  mesh  size)  was  conducted  for  two  sphere  problems. 
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TABLE  3.  RELATIVE  ERROR  (PERCENT)  AND  CP  TIME  (SECONDS) 
FOR  EXAMPLE  12  FROM  REFERENCE  2 


BANDWIDTH 

ITERATIONS 

6 

10 

14 

18 

22 

1 

*t7.09 

22.36 

10.70 

12.1*7 

2 

23.98 

*♦.95 

1.25 

1.67 

HI 

3 

12.**9 

1.17 

.22 

.23 

HI 

l* 

6.53 

.28 

5 

3.1*2 

91 

6 

1.79 

1 

7 

.9** 

TIME 

0.738 

0.477 

0.1*53 

0.577 

0.591 

EFF  (CP  TIME) 

5.7 

8.8 

9.2 

7.2 

7.1 

EFF  (MO'S) 

5.3 

8.5 

8.8 

6.8 

6.5 

L 


The  reader  is  referred  to  reference  k for  the  details  and  to  reference  8 
for  a summary.  The  resulting  guidelines  call  for  segment  lengths  somewhat 
less  than  0.2A  as  a maximum  and  a wire  radius  of  about  0.025A.  Predicted 
•.  far  field  patterns  and  bistatic  scattering  cross  sections  for  models 

constructed  with  these  guidelines  were  consistently  in  good  agreement 
1 with  the  independent  results  from  the  literature. 

If 

The  modeled  objects  in  this  study  all  had  rotational  symmetry, 
i The  models  were  constructed  with  many-fold  rotational  symmetry.  Both 

1 this  symmetry  and  the  iterative  method  were  used  in  the  solution  process. 

Due  to  symmetry  operations  performed  on  the  matrix,  the  iterative  method 
was  not  highly  efficient  in  this  study  because  large  submatrices  and  band- 
' widths  were  required.  This  is  not  especially  important,  because  the 

only  significant  cost  for  such  problems  is  in  generating  part  of  the 
impedance  matrix. 

Reference  5 includes  a number  of  example  problems  with  N 
f ranging  up  to  1000.  The  solutions  were  mostly  obtained  without  using 

symmetries.  Banded  matrix  iteration  is  shown  to  be  a useful  solution 
method  for  large  problems. 

Most  of  these  example  problems  were  chosen  because  the  predicted 
t current  distribution  or  far  fields  could  be  compared  to  independent 

results  from  the  literature.  Few  problems  of  the  desired  electrical 
size  range  have  been  solved  for  current  distributions,  so  most  of  the 
comparisons  are  to  theoretical  far  fields.  As  noted  previously,  reasonably 
good  far  fields  can  be  obtained  from  inaccurate  currents.  Comparison  to 
theoretical  far  fields  is  not  the  best  validation  for  choice  of  model 
parameters . 
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Comparisons  of  shapes  of  grid  current  distributions  and  sur- 
face current  distributions  are  possible  in  some  cases.  The  wires 
usually  must  be  oriented  on  the  surface  along  the  natural  coordinates, 
because  surface  currents  are  computed  by  component  along  those  coordi- 
nates. A rectangular  plate  modeled  by  a square  mesh  grid  is  an  obvious 
example.  For  a sphere,  the  natural  wire  orientations  are  along  lines  of 

• latitude  and  longitude.  For  plane  wave  scattering,  surface  currents  are 
known  for  the  E-  and  H-planes.  If  the  plane  wave  is  incident  along  the 
polar  axis,  H-plane  currents  are  ^-directed  and  can  be  compared  to 
currents  on  wires  crossing  the  H-plane  and  oriented  along  lines  of 
latitude.  E-plane  surface  currents  cannot  be  compared  to  currents  on  a 
set  of  wires  along  a line  of  longitude,  because  these  currents  break  up 
and  run  through  the  pole  on  a variety  of  wire  paths.  If  the  plane  wave 
is  incident  along  a tine  normal  to  the  polar  axis  and  the  electric  field 

f is  in  the  equatorial  plane,  both  comparisons  are  possible.  In  the  E- 

plane,  the  grid  is  essentially  rectangular.  Currents  are  ^"directed  in 
the  H-plane,  and  can  be  compared  to  currents  on  wires  along  lines  of 
latitude  except  at  the  poles.  Such  a comparison  is  shown  in  figure  14 

• for  a sphere  with  ka  equal  to  4.7  (example  13  from  reference  5).  The 
model  used  996  wire  segments  for  a surface  area  of  7 square  wave- 
lengths. The  predicted  bistatic  scattering  cross  section  was  in 
excellent  agreement  with  the  exact  theoretical  results  of  King  and 

Wu  preference  14). 

In  general,  the  iterative  process  yields  good  efficiencies  for 
planar  and  near-planar  wire  grids  when  the  strip  numbering  scheme  is 
used  and  the  bandwidth  is  chosen  from  figure  9.  It  also  works  well 
tor  cylinders  with  diameters  less  than  about  a half-wavelength.  For  fat 
. vlinders  and  spheres,  larger  bandwidths  are  required  for  convergence. 

. v '.tths  chosen  from  figure  9 for  these  problems  yield  pseudocon- 

e or  divergence.  No  rule  for  bandwidth  selection  in  such  cases 
ice  the  study  of  such  problems  was  brie".  It  is  possible 
’jl  modes  contribute  to  this  problem;  example  15  of  reference 
at  to  at  to  investigate  this  possibility. 
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FRONT  ANGLE  (deg)  REAR 


Figure  14.  Amplitude  of  Current  on  a Sphere  for  Plane  Wave  Scattering, 
ka  = 4.7.  Example  13  of  Reference  5 (Surface  Currents  K 
from  Reference  14;  Solid  Line  for  E-Plane,  Dashed  for 
H-Plane.  Wire  Grid  Currents  I;  © for  E-Plane,  □ for  Near 
H-Pl ane. ) 
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As  an  example  of  the  application  of  the  iterative  process  in 
larger  problems,  example  6 of  reference  5 is  reproduced  here. 

EXAMPLE  6 OF  REFERENCE  5 

A flat  square  plate  is  modeled  as  shown  in  table  A.  A plane 
wave  with  normal  incidence  is  scattered  by  the  wire  grid.  The  incident 
E-field  is  parallel  to  one  edge  of  the  plate. 

In  example  6A,  the  wire  segments  are  numbered  in  strips  at 

right  angles  to  the  incident  E-field.  The  plate  edge  length  is  2X. 

Figure  9 indicates  a bandwidth  for  peak  efficiency  corresponding  to 

a distance  of  O.jA  for  an  object  of  length  2A.  With  segment  lengths  of 

0.125A,  a bandwidth  of  132  includes  interactions  to  a minimum  distance 

of  0.5A.  Table  5 shows  the  convergence  data.  Convergence  is  so  rapid 

that  a smaller  bandwidth  would  probably  provide  better  efficiency.  The 

normalized  bistatic  cross  sections  in  the  E-  and  H-planes  are  shown  in 

2 

figures  15  and  16.  The  backscatter  cross  section  o/X  is  220.  The 
pivot  ratio  was  20.0.  The  matrix  fill  time  was  608  seconds  and  the 
solution  time  was  138  seconds. 

Example  6B  is  identical  to  example  6A  except  for  orientation 
of  the  incident  field,  wrich  is  parallel  to  the  X-axis.  This  is  equiva- 
lent to  retaining  the  orientation  parallel  to  the  Y-axis  and  numbering 
in  strips  normal  to  the  X-axis.  The  two  are  physically  equivalent 
problems.  The  convergence  data  are  shown  in  table  6.  The  solutions 
for  the  two  examples  were  the  same  within  about  1 percent  in  the  regions 
of  large  current. 

Currents  parallel  to  the  field  are  shown  in  figures  17(a),  17(b), 
and  17(c)  for  wires  along  the  edge,  next  to  the  edge,  and  in  the 
interior  of  the  grid.  The  edge  wires  carry  a larger  current  as  expected. 
The  wires  next  to  the  edge  carry  slightly  less  current  than  interior 
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DESCRIPTION  OF  EXAMPLE  6 FROM  REFERENCE  5 


ANGLE  FROM  INCIDENCE  (deg) 


Figure  15.  Normalized  E-Plane  Bistatic  Cross  Section 
for  a Wire  Gridded  Plate,  Example  6A 


ANGLE  FROM  INCIDENCE  (deg) 


Figure  16.  Normalized  H-Plane  Bistatic  Cross  Section 
for  a Wire  Gridded  Plate,  Example  6A 
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Currents  normal  to  the  E-field  are  zero  along  the  centerline 
wires  due  to  symmetry.  They  are  small  everywhere  compared  to  currents 
parallel  to  the  field,  and  are  largest  at  the  edges.  Figure  18  shows 
the  approximate  edge  current  normal  to  the  field.  These  small  currents 
are  subject  to  considerable  error  due  to  the  1 percent  convergence 
criterion.  The  upturn  at  the  corners  is  due  to  the  difference  in  bound- 
ary conditions  between  plates  and  wire  grids. 

Comparison  of  tables  5 and  6 shows  that  the  error  in  Xj  is 
considerably  larger  for  example  6B,  although  the  final  efficiencies  are 
not  much  different.  The  higher  initial  error  is  due  to  the  direction  of 
the  E-field  (and  hence  the  currents)  relative  to  the  orientation  of  the 
strips  used  in  segment  numbering.  Interactions  between  segments  at 
moderate  distances  are  largest  when  the  segments  are  parallel  to  each 
other  and  normal  to  the  line  between  their  centers.  In  example  6A,  the 
banded  matrix  includes  distant  interactions  between  parallel  segments 
with  large  currents  and  between  col  linear  or  nearly  col  linear  segments 
with  essentially  no  current.  In  example  6B,  the  opposite  is  true. 

Hence,  although  exactly  the  same  matrix  operations  are  involved  in  the 
two  cases,  the  excitation  affects  the  results  considerably. 

Example  6C  differs  from  example  6A  only  in  the  wavelength. 

With  A equal  to  1.176,  the  quantity  kW/2  or  ttW/A  is  equal  to  the  square 
root  of  28.  Ufimtsev  (reference  15)  has  investigated  plane  wave  scattering 
from  an  infinite  strip  with  this  parameter  value  tor  the  strip  width.  For 
normal  incidence  and  E across  the  strip  (H-polarizat ion) , he  gives  the 
normalized  far  field  amplitude  shown  in  figure  19.  For  E parallel  to 
the  strip  (E-polarization),  he  gives  the  no  malized  far  field  amplitude 
shown  in  figure  20  by  a solid  line.  With  the  same  normalization,  the 
far  field  in  the  E-plane  for  the  wire  grid  is  identical  to  figure  19 
within  graph  reading  limitations.  The  far  field  in  the  H-plane  differed 
from  the  strip  results  as  shown  by  dots  in  figure  20. 
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DISTANCE  ACROSS  PLATE 


Figure  18 


Wire  Currents  Along  Edge  Normal  to 
The  E-Field,  Example  6 


Figure  19.  Far  Field  for  Plane  Wave  Scattering  from  a 
Strip,  H-Pol ar i zat ion , Reference  15 


Figure  20.  Far  Field  for  Plane  Wave  Scattering  from  a 
Strip,  E-Polarization  (Line  from  Reference 
15,  Points  from  Wire  Grid) 


81 


Convergence  data  for  example  6C  are  shown  in  table  7.  The 
pivot  ratio  was  9-7*  The  matrix  fill  time  was  592  seconds,  and  the 
solution  time  was  136  seconds. 


E.  RECOMMENDATIONS  FOR  FUTURE  EFFORTS 


The  development  of  an  EM  analysis  code  using  the  GEMACS  architec- 
ture should  establish  a methodology  for  future  efforts. 

The  present  GEMACS  code  is,  in  essence,  a data  base  manipulation  tool 
which  also  supports  certain  types  of  operations  on  the  data  involved. 

The  current  operations  result  in  the  EM  analysis  of  those  structures  and 
antennas  which  can  be  modeled  using  the  thin-wire  approximation  to  the 
EFIE  with  a sine  + cosine  + pulse  expansion  function.  Additional  capa- 
bility in  the  method  of  moments  areas  should  be  added  which  would  allow  a 
consistent  analysis  of  the  modeling  requirements  based  on  different 
expansion  and  weighing  functions. 

Additionally,  incorporation  of  the  MFIE  (Magnetic  Field  Integral 
Equation)  would  extend  the  capability  considerably  and  allow  investigation 
of  the  BMI  technique  with  a different  type  of  matrix. 

Besides  expansion  of  the  method  of  moments  oriented  solution  tech- 
nique, diffraction  theory  methods  could  be  included.  If  properly  done, 
this  would  support  analysis  using  method  of  moments  derived  current 
distributions  with  diffraction  theory  field  computations. 

With  additional  software  for  matrix  manipulation,  the  diakoptic 
solution  process  could  be  investigated  to  determine  its  applicability 
to  field  coupled  problems.  This  technique  is  already  applied  to  distributed 
power  networks  and  mechanical  systems  with  a great  deal  of  success. 

Besides  increasing  the  capability  of  GEMACS  to  solve  a larger  class 
of  problems,  there  are  some  software  improvements  which  could  greatly 
reduce  the  size  and  run  time  of  the  program.  These  improvements  would 
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not  correspond  to  ANSI  FORTRAN  IV;  however,  the  resultant  code  would 
still  run  on  a large  number  of  second  and  third  generation  computers. 

The  two  most  obvious  improvements  are  random  access  I/O  and  dynamic  core 
allocation.  Both  of  these  operations  require  software  peculiar  to  the 
host  machine  and  are  generally  available. 
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APPENDIX  A 

ANTENNA  SYSTEM  EFFECTIVENESS  AS  AFFECTED  BY  ENVIRONMENT 

The  effectiveness  of  an  antenna  system  can  be  adversely  affected  by 
its  environment.  As  one  example,  consider  the  frequently  encountered 
problem  in  antenna  siting  of  realistically  including  the  ground  effects 
upon  the  antenna  behavior.  The  antenna  properties  can  be  influenced  by 
the  ground  in  essentially  two  ways:  (1)  by  modifying  the  antenna  current 

distribution  through  the  near-field  interaction;  (2)  and  by  changing  the 
far-field  radiation  pattern  from  its  free-space  value.  Even  in  the  case 
of  a perfectly  conducting  ground,  where  image  theory  is  applicable  in  a 
straightforward  way,  the  mutual  interaction  between  the  antenna  and  its 
ground  image  must  be  allowed  for  if  the  coupling  is  significant. 

The  antenna  engineer  is  oftentimes  faced  with  the  task  of  designing 
a system  with  certain  character i s t i cs  in  a particular  ground  environment. 
The  rigorous  approach  presented  by  Sommerfeld  (1964)  is  somewhat  cumber- 
some and  lengthy  even  for  computer  evaluation.  An  accurate  but  approxi- 
mate method  which  surmounts  the  numerical  difficulties  associated  with 
the  Sommerfeld  approach  is  used  in  the  GEMACS  program.  A tested  integral 
equation  approach  for  free-space  antenna  analysis  was  modified  in  a 
straightforward*  manner  to  provide  an  accurate  and  efficient  numerical 
procedure  for  handling  the  antenna/ha  1 f -space  problem.  Fresnel  plane- 

wave  reflection  coefficients  are  used  to  allow  modeling  of  the  ground. 

C 

A.  THE  DIFFERENTIAL  FIELD  OVER  GROUND 

1 . Vertical  Source  Segment 

The  Hertzian  potential,  dt^,  of  a vertically  directed  electric 
current  filament  of  strength  Id  located  in  free  space  above  a half- 
space of  permittivity  and  conductivity  defined  by  z <_  0 is  given 
by  (Sommerfeld,  1964) 

dnv(r)=  zdir^(r)  = z P I ( r')  d 8,{  g^  + g.  - g^}  (Eq.  A- 1 ) 
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2. 


Horizontal  Source  Segment 


The  Hertzian  potential  dir^  of  a horizontal  differential  electric 
current  source  ld£  pointing  at  an  angle  a with  respect  to  the  x axis  at 
coordinates  x' , y' , z'  above  a half-space  defined  by  z £ 0 of  relative 
permittivity  e is  given  by 


(Eg.  A-2) 


dirh  (r)  = Pl(r')d*  [9d  - g.  + gssJ  s 


+ g z = dir  , x + dir  , y + drc  , z 
3sz  xh  yh  1 zh 


where  g^  and  g.  are  given  by  eguation  A-l  and 

CO 

g =2  /*  J Up")  e-Az+Z']  4-dX 

ss  o P+u,. 


9sz  = i?  ccs 


f J'  Up")  e“u(z+z")  -■  ^ ■ ■ X2  dX 

o eEp+pE 


p"  = [(x-xu  + (y-y')2] 


s = x cos  a + y sin  a, 


as  depicted  in  figures  A-l  and  A-2,  and  is  the  derivative  of  Jg  with 
respect  to  its  argument.  It  should  be  observed  that  the  horizontal 
antenna  produces  not  only  a Hertzian  potential  aligned  with  the  current 
filament,  but  a vertical  component  as  well  as  a result  of  the  half-space 
interaction.  Note  also  the  change  in  sign  of  the  image  term  g.  compared 
with  the  vertical  source. 
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For  both  of  the  above  cases,  an  expression  for  the  differential 
electric  field,  de  of  the  source  is  given  by 

de  (7)  = k2  dir  (7)  + VV  • dn  (7)  (Eq.  A-3) 

where  V and  V • operate  on  the  unprimed  or  observation  coordinates,  and 
for  simplicity,  the  subscripts  v and  h which  denote  the  source  orienta- 
tion are  omitted.  Note  that  three  terms  in  the  bracket  of  equation  A-l 
correspond  respectively  to  the  direct  radiation  field,  the  perfectly 
conducting  half-space  image  field,  and  the  Sommerfeld  integral  correction 
for  finite  ground  conductivity.  The  latter  two  terms  account  for  the 
ground  or  interface-reflected  contribution  to  the  total  field. 

B.  INTEGRAL  EQUATION 


In  order  to  determine  the  electric  field  from  a source  consisting 
of  an  arbitrary  collection  of  wire  segments,  it  is  necessary  to  integrate 
these  differential  fields  over  the  source  structure 


E 


dir  (?)  + VV 


dir  (7) 


(Eq.  A--4) 


In  order  to  determine  the  current  on  a wire  antenna  located 
over  a ground  plane,  the  following  integral  equation  must  be  solved 


J 

sr.  ’ E (rJ  " s„  * E (O  = 0 

o o o o 


(Eq.  A-5) 


where  s and  r are  defined  in  subsections  A. 1 and  A. 2. 
o o 

This  integral  equation  is  rigorously  correct  (within  the  restrictions 
imposed  by  the  thin-wire  approximation)  for  a wire-antenna  located  in 
free-space  over  a ground  of  arbitrary  permittivity  and  conductivity. 

The  resulting  integral  equation  thus  has  a portion  which  involves  the 
double  integration  over  and  X,  a procedure  that  in  general  cannot  be 
analytically  carried  out. 


A-6 


In  view  of  this,  it  is  advantageous  to  make  use  of  reasonable 
approximations  to  provide  a method  which  is  simpler  to  use  and  which 
does  not  significantly  degrade  the  numerical  accuracy  of  our  results. 

C.  THE  REFLECTION  COEFFICIENT  APPROXIMATION 

One  approximate  solution  for  the  Sommerfeld  integral  employs 
asymptotic  expansions  for  the  integral  which  are  valid  over  various 
ranges  of  the  parameter  values.  This  particular  method  makes  use  of  the 
method  of  stationary  phase  and  double  saddle-point  integration.  The 
results  from  using  this  approach  include  closed  form  expressions  for 
the  fields  of  an  elementary  source,  electric  or  magnetic,  useful  over 
various  ranges  of  observation  point  location  relative  to  the  source, 
and  involving  the  electrical  parameters  of  the  two  media. 

A comprehensive  survey  of  this  general  methodology  is  given  by 
Banos  (1966),  and  the  numerical  accuracy  of  these  approximations  is 
examined  by  Siegel  (1970).  Siegel's  work  is  especially  useful  in  that 
he  is  able  to  establish  ranges  of  observation  distance  for  which  the 
various  expressions  are  valid,  but  unfortunately  also  shows  that  gaps 
exist  in  general  over  which  the  approximations  do  not  provide  acceptable 
accuracy.  This  is  true  in  particular  for  observation  points  located 
close  to  the  source,  so  that  the  integral  equation  A-3,  which  requires 
this  close  proximity  of  observation  and  source  points  is  not  amenable  to 
using  these  approximate  expressions. 

Other  formulations  of  this  general  problem,  suitable  for  possibly 
more  physically  based  approximations,  have  been  investigated.  The 
reader  is  referred  to  the  work  by  Feynberg  (1967),  where  a surface 
integration  over  the  interface  induced  sources  is  discussed,  with  many 
other  aspects  of  the  problem,  and  the  conference  proceedings  edited  by 
Wait  (1969a).  The  latter  provides  a number  of  interesting  treatments, 
including  one  based  on  the  compensation  theorem  (Surtees,  1969)  and  a 
comparison  of  several  formulations  for  antenna  impedance  (Wait,  1969b). 
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In  developing  an  approximate  solution  to  equation  A-5,  we  chose  to 
proceed  on  a primarily  physical  basis,  and  further  to  seek  an  approach 
which  is  well  suited  to  the  integral  equation  treatment  developed  for 
I the  free-space  antenna  problem.  We  have  noted  that  the  first  two  terms 

* in  the  Green's  function  of  equations  A-l  and  A-2  are  the  direct  or  free- 

space  radiation  and  the  perfectly  conducting  half-space  or  image  radia- 
tion, respectively.  The  image  term  represents  the  contribution  of  a 
specularly  reflected  ray  from  a perfectly  conducting  half-space,  the 
electromagnetic  equivalent  of  an  optically  perfect  mirror  in  optics. 

This  ray  optics  specular  component  is  thus  rigorously  correct  for  the 
perfect  conductor  for  which  the  Sommerfeld  correction,  gs>  becomes  a 
zero.  Further,  the  X-integration,  as  shown  in  equations  A-l  and  A-2  can 
be  analytically  performed  for  the  image  term  (and  the  direct  term  as 
well)  so  that  there  is  no  significant  increase  in  difficulty  in  treating 
an  antenna  located  over  a perfectly  conducting  half-space  over  analyzing 
the  same  antenna  in  free-space. 

The  advantage  of  extending  the  ray  optics  approach  to  the  finitely 
conducting  half-space  is  immediately  evident.  If  a specularly  reflected 
ray  can  be  used  in  this  case  as  well  to  account  for  the  entire  ground- 
reflected  contribution  (g.  and  g^  terms)  to  the  field,  then  it  may  be 
possible  to  largely  circumvent  the  numerical  difficulties  presented  by 
the  Sommerfeld  integral.  Two  problems  arise  in  accomplishing  this  --  one 
is  the  determination  of  the  appropriate  reflection  coefficient  for  the 
half-space,  and  the  other  is  that  of  specifying  the  specular  ray  contri- 
bution itself.  These  problems  are  both  brought  about  because  the  fields 
involved  are  non-plane  wave  and  the  observation  point  can  in  general  lie 
in  the  near-field  of  the  source  where  the  higher  order  terms  may  be  of 
importance. 

Some  insight  may  be  gained  by  viewing  the  interface  as  a surface 
distribution  of  induced  sources,  as  discussed  by  Feynberg.  The  ground- 
reflected  contribution  to  the  total  field  produced  by  a given  infinitesi- 
mal source  is  then  obtained  from  a surface  integration  over  these  induced 
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sources.  Practically  speaking,  this  surface  integration  can  be  truncated 
to  include  a region  of  only  limited  extent.  Upon  applying  the  principle 
of  stationary  phase,  this  truncated  region  (or  essential  region  as 
denoted  by  Feynberg)  can  be  shown  to  be  of  generally  elliptical  shape 
and  centered  in  some  sense  about  the  specular  reflection  point.  If  the 
area  of  this  essential  region  is  allowed  in  the  limit  to  shrink  to  zero 
centered  on  the  specular  point,  only  the  specular  ray  itself  will  be 
obtained.  Thus,  a ray  optics  approach  based  on  a specular  ray  amounts 
essentially  to  collapsing  the  entire  interface  induced  source  distribu- 
tion to  a single  point  from  which  the  ground  reflected  field  at  a partic- 
ular observation  point  to  a given  infinitesimal  source  appears  to  originate. 

This  is  equivalent  to  the  use  of  an  image  source  to  account  for  the 
reflected  field.  The  principle  applies  to  extended  sources  as  well,  by 
integrating  over  the  corresponding  specular  point  distribution  of  equiva- 
lent image  sources.  Since  this  picture  of  the  reflection  process  is  of 
course,  exact  for  the  perfectly  conducting  half-space,  its  extension  to 
the  finitely  conducting  problem  with  a minimum  of  modification  appears 
expedient.  To  obtain  the  functional  dependence  of  the  total  ground 
reflected  field  upon  the  geometrical  parameters,  it  follows  that 
equation  A-l  with  the  g.  term  only  should  be  used.  Complete  specular 
reflection  of  the  field  components  of  all  orders  is  thus  presumed  (unit 
magnitude  reflection  coefficients  assumed).  The  effect  of  finite  ground 
conductivity  consequently  enters  below  only  via  a modified  reflection 
coefficient  which  multiplies  the  image  fields,  and  thus  effectively 
alters  the  image  strength. 

Again,  the  simplest  extension  beyond  the  perfectly  conducting  case 
which  will  introduce  the  ground  electrical  parameters  in  an  apparently 
meaningful  way  is  to  be  preferred.  Certainly  the  simplest  reflection 
coefficient  which  could  be  employed  would  be  one  based  on  plane  waves, 
the  so-called  Fresnel  reflection  coefficients.  Since  these  are,  however, 
polarization  dependent,  and  the  vertical  and  horizontal  field  components 
reflect  in  different  ways  as  well,  the  ground  reflected  field  must  be 
decomposed  into  the  usual  TE  (transverse  electric)  and  TM  (transverse 
magnet  i c modes ) . 
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This  is  readily  accomplished  by  breaking  the  image  field,  or  g. 
terms  of  equation  A-**,  into  electric  field  components  lying  in  TM  and 
perpendicular  to  TE  the  vertical  plane  containing  the  source  and 
observation  points  (incidence  plane).  The  TM  components  are  further 
identified  as  a horizontal  or  vertical  relative  to  the  horizontal  inter- 
face. Each  of  these  field  components  is  then  multiplied  by  the  corres- 
ponding plane  wave  electric  field  reflection  coefficient,  using  as  the 
angle  of  incidence  the  angle  made  by  the  specular  ray  connecting  the 
source  and  observation  points  (midpoints  of  the  source  and  observation 
segments,  respectively).  This  procedure  thus  neglects  the  possible 
angle  of  incidence  variation  resulting  from  a non-zero  length  source 
segment  for  which  the  actual  specular  point  distribution  would  be  the 
path  on  the  interface  traced  by  a straight  line  drawn  from  the  observa- 
tion point  and  scanned  over  the  image  segment.  Finally,  the  total 
tangential  observation  point  field  is  obtained  as  the  sum  of  the  direct 
and  ground  reflected  fields. 

The  approach  described  above  is  a straightforward  extension  to  the 
E F I E 1 for  free-space,  or  to  the  perfectly  conducting  half-space.  Yet, 
as  will  be  demonstrated  by  the  numerical  results  which  follow,  it  provides 
surprisingly  good  results  for  the  analysis  of  wire  antennas,  over  a 
half-space  with  widely  varying  electrical  parameters.  Since  the  method 
is  based  on  a surface  reflection  coefficient  at  the  specular  point,  it 
is  readily  adapted  to  a horizontally  stratified  half-space,  and  to  a 
half-space  with  slowly  varying  properties  along  the  interface. 

For  reasons  of  simplicity,  it  is  assumed  that  the  antenna  is  located 
in  free-space  above  a horizontal  ground  with  the  x-y  plane  forming  the 
interface.  Since  any  wire  segment  can  be  specified  in  terms  of  three 
orthogonal  components,  it  is  only  necessary  to  determine  the  differential 
fields  due  to  a vertical  Hertzian  source  and  an  arbitrarily  oriented 
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horizontal  source.  The  following  development  is  limited  to  the  vertical 
source.  Subsequent  sections  will  deal  with  the  arbitrarily  oriented 
horizontal  wire  segments  in  a similar  fashion. 

The  reflection  coefficient  approximation,  the  name  by  which  we 
characterize  the  modified  image  method  outlined  above,  represents  only 
a simple  extension  to  the  EFIE  for  the  case  of  a wire  antenna  located 
over  a perfectly  conducting  half  space.  To  provide  some  theoretical 
justification  for  the  reflection  coefficient  approximation  it  will  be 
compared  with  the  rigorously  correct  Sommerfeld  integral  formulation. 

1 . Vertical  Filament 

It  will  be  convenient  in  comparing  the  Sommerfeld  and  reflection 
coefficient  approaches  (and  no  loss  of  generality  results)  to  deal  with 
the  differential  fields  dE  which  result  from  an  incremental  current 
filament  dl  located  on  the  z axis.  In  addition,  the  principal  polarized 
fields  are  simply  the  cylindrical  coordinate  system  components.  Let 


then,  from  equation  A-l  and  A-3 

d E (p,  z)  = Q [z  + v-y  (gd  + g.  - gs) 

The  totai  differential  field  can  be  written  as  the  sum  of  dE  the 

. ^ 

direct  component  (due  to  g^)  and  dE  , the  reflected  component  (due  to 
g.  - gs)  so  that 

dl  = dl  D + dl  R 

As  noted  above,  the  reflection  coefficient  approximation  involves 

R 

computing  dE  in  terms  of  the  perfect  ground  contribution;  let  us 

denote  the  latter  image  component  by  d E (due  to  g.).  Then  dE  is 

- R 1 

approximated  by  dE  where 


dE  R 2J  dE  R = r dE  1 
v 

where  r is  the  Fresnel  reflection  coefficient  matrix  for  the  vertical 
v 

source,  given  by 


0 


= R"  pp  + R"  zz 
I V 


0 R"  \ 

v J 


(Eq.  A-6) 


Note  that  r is  a geometrical  function  of  the  relative  locations  of  the 
source  and  observation  points  only.  The  superscript  denotes  the  electric 
field  polarization  relative  to  the  incidence  plane  (in  this  case  the  p, 
z plane)  and  the  subscript  denotes  the  field  component  relative  to  the 
interface.  While  only  the  z-  or  vertical  field  is  required  for  the 
integral  equation  solution,  the  horizontal  component  is  required  for  the 
far  field  calculation,  and  is  included  for  completeness.  Of  course,  the 
differential  electric  fields  must  be  similarly  decomposed  as 

II  II 

R R,  * , r R, 

dE  = dEu  ’ p + d E ’ z 

H v 


with  corresponding  expressions  for  the  other  fields.  Note  that  the 

vertical  dipole  produces  only  z-  and  p-  directed  electric  fields  since 

the  cf>-  components  of  the  various  Vg  terms  are  zero. 

The  coefficients  in  T are  the  standard  Fresnel  plane  wave  reflection 

v 

coefficients  given  by 


e cos  0 - - sin2  0 ,, 

R"  = — = - R ,, 

V „ / ; — t — — H 

cos  0 + ✓ - sm^  0 


where  the  angle  of  incidence  relative  to  the  interface  normal  is 
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Observe  that  the  plane-wave  reflected  electric  fields  are  related  to 
the  incident  field  in  the  same  manner  as  equation  A-6,  with  ||  corres- 
ponding to  the  TM  wave  polarization.  Note  that  the  reflection  coeffic- 
ients are  constant  for  the  self-fields  of  a single  vertical  antenna, 
where  the  angle  of  incidence  is  effectively  zero,  since  x = x'  and 
y = y'  everywhere  on  the  antenna.  For  radially  separated  vertical  anten- 
nas, this  equality  no  longer  holds,  however,  and  the  variation  of  the 
matrix  with  observation  and  source  coordinates  must  be  taken  into  account. 

2.  Arbitrarily  Oriented  Horizontal  Filament 

For  the  case  of  an  arbitrarily  oriented  horizontal  segment,  the 
reflection  coefficient  approximation  requires  that  the  electric  field  at 
a given  observation  point  due  to  a differential  source  (or  a single  antenna 
segment)  be  decomposed  into  components  which  correspond  to  the  transverse 
magnetic  (TM)  and  transverse  electric  (TE)  polarizations  of  a plane  wave 
relative  to  the  incidence  plane.  It  is  thus  necessary  in  general  to  find 
the  three  orthogonal  components  of  the  field,  two  of  which  lie  in  the 
plane  of  incidence  (for  a non-zero  length  segment,  the  source  point  is 
taken  to  be  the  segment  center),  and  one  which  is  orthogonal  to  it.  The 
two  in-plane  field  components  are  the  radially  and  vertically  (or  z) 
directed  parts  of  the  field,  while  the  orthogonal  component  is  the  azi- 
muthally  directed  part,  relative  to  a cylindrical  coordinate  system  on 
whose  vertical  or  z-axis  the  source  is  located.  Since  the  total  field 
is  obtained  as  the  summation  of  differential  fields  due  to  a distribution 
of  sources  which  will  not  in  general  have  a common  incidence  plane,  it  is 
necessary  to  perform  the  field  integration  in  a convenient,  absolute 
coordinate  system,  such  as  cartesian,  or  cylindrical,  with  respect  to 
which  the  source  distribution  is  specified. 
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It  is  consequently  expedient  to  express  the  differential 
field  directions  relative  to  the  source-incidence  plane,  and  to  then 
transform  this  field  into  components  relative  to  the  absolute  reference 
system.  Let  the  in-plane  unit  vectors  which  are  horizontal  and  vertical 
relative  to  the  interface  be  s|(  and  z respectively  and  let  the  orthogonal 
unit  vector  to  the  plane  be  s^  . The  three  field  components  at  observa- 
tion point  x,y,z  (or  p,  4>,  z)  due  to  a source  at  x' , y',  z'  (or  p',  <)>',  z') 
can  be  expressed  as 


dE  = de  x + de  y + de  z = de 
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Note  that  <(>"  is  the  angular  coordinate  of  the  observation  point  relative 
to  the  source,  and  y is  the  angular  separation  between  the  incidence 
plane  and  a radial  vector  through  the  observation  point.  For  the  special 

case  of  the  vertical  source  de  is  zero. 

x 

As  a final  point  concerning  the  reference  directions  for  the 
electric  field  calculation,  it  is  pertinent  to  mention  that  the  direct, 
image  and  Sommerfeld  terms  of  the  Green's  function  produce  only  longi- 
tudinal and  radial  differential  fields  relative  to  the  source  direction. 
This  is  due  to  the  azimuthal  independence  of  a filamentary  source  and 
its  resultant  fields.  Thus  the  differential  field  is  most  conveniently 
found  for  an  arbitrarily  directed  source  in  the  source-oriented  longi- 
tudinal and  radial  directions,  and  then  subsequently  transformed  into 
the  appropriate  components  required  for  the  integral  equation  calcula- 
tion. For  computer  implementation  it  is  generally  most  convenient  to 
arrange  the  overall  numerical  procedure  entirely  in  terms  of  cartesian 
components.  In  the  following  discussion,  it  will  be  more  illuminating 
in  some  instances  to  use  the  sM,  s and  z directions  in  demonstrating  the 
application  of  the  reflection  coefficient  approximation. 

Let  the  radial  and  longitudinal  electric  fields  referred  to 
the  source  (or  image)  oriented  cylindrical  coordinate  system  (p,  6,  Z,) 

be  denoted  by  dE  and  dE„  respectively.  These  fields  are  obtainable 

p l 

f rom 


dE  = dE  p + dE  l 
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(Eq.  A-8) 


dE  = dE  p + dE  l 
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for  a segment  of  length  A where  g is  the  Green's  function  e /r  with  r 
the  distance  from  the  observation  point  to  the  source  (or  image).  In  the 
source  centered  system 

r = /pz  + l 


with  p,  6,  and  l the  observation  point  coordinates.  As  previously  noted, 

it  is  possible  to  analytically  evaluate  the  dE  and  dE  terms  of  equation 

P 

A-8  for  a cosine  or  sine  current  variation,  as  well  as  the  dE^  component 
for  a current  independent  of  l'. 

We  can  now  transform  from  the  p and  l field  components  to  the 
cartesian  components  using 


x 

y 


sin  a sin  6 
-cos  a sin  6 

cos  6 


sin  a cos  6 
-cos  a cos  6 

-sin  6 


cos  a 
s in  a 


where  6 is  measured  from  the  vertical.  Continuing,  we  then  proceed 
to  the  incidence-plane  image  field  components  which  are  expressible  in 
terms  of  the  unit  vector  sM  and  sA  via  the  transformation  given  by 
equation  A~7.  The  required  image  field  components  (which  we  denote  with 
a superscript  I and  which  come  from  equation  A-8  with  g = g.  relative 
to  the  incidence  plane)  can  be  put  in  the  form 


) 
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f 


m 

” i 

cos  <t>‘ ' sin  0 

dE1! 

= 

-sin  (ji"  cos  <f>"  0 

dE1 

0 0 1 

z_ 

sin  a sin  6 sin  a cos  <$  cos  a 
-cos  a sin  6 -cos  a cos  6 sin  a 
cos  6 -sin  6 0 


dE 


dE 


dE 

L Z. 


where 


= tan 


(&) 

[(x-x')  sin  n + (y-yOcoscj 

J 

I 1 + cot  a tan  <j>"  ~1 

(_  sin  a (x-x""  ) (z-z')  J 

p = [(x-x^)2  + (y-y')2  + (z-z')2  - l2  J 

l = (x-x")  cos  a + (y-y')  sin  a 


6 = tan 


= tan 


Note  that  the  direct  field  components  and  the  image  components  as  obtained 
from  equation  A-8  differ  only  because  z,'  is  the  positive  for  the  former, 
and  negative  for  the  latter. 

Since  the  dE^  component  of  electric  field  is  zero,  the  above 
simplifies  to 


dE'„  - 

s i n 

6 

s i n 

( ot — 4>" ) 

< - 

-s  i n 

6 

cos 

(a-4>") 

Cl 

m 

N — 

II 

cos 

6 

Q. 

m 

D — 

Upon  mu  1 1 i pi i cat ion_of  the  e/,,  E ^ and  E^  components  by  the  reflection 
coefficient  matrix  f,  where 
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l 


where 


s„  s„  + Rh  s s + Ry 


z z 


cos  0 - /e_  - sin^  9 

o E 


cos  0 + /er  - sin'^  0 

o E 


we  obtain  the  desired  reflected  field  components,  which  we  denote  by 
superscript  R, 


dE,, 


" rh  [sin 

J = RH[Sin 


<5  sin  (a-<j>")  dE*  + cos  (a-<f>“)  dE! 

P ^ 


6 cos  (a-<t>")  dE^  + sin  (a-4>")  dE^ 


Dll  I 

dE  = R cos  6 dE1 
z v p 


As  the  final  step  in  the  procedure  of  approximating  equation 
A-5,  add  the  direct  and  reflected  fields  (transformed  back  to  cartesian 
components)  due  to  the  given  source  segment,  evaluated  at  an  observation 
point  at  another  ar.'er’a  seg-.ent,  and  take  the  dot  product  with  the  obser 
vat  ion  segment  tangent  vector. 

(1)  Pass  a vertical  plane  through  the  center  of  the  source 
segment  and  the  observation  point.  (Plane  of  Incidence) 

(2)  Determine  the  angle  of  incidence  at  the  specular  point. 

(3)  Decompose  the  electric  field  into  components  parallel 
and  perpendicular  to  the  plane  of  incidence. 

(4)  Apply  the  appropriate  reflection  coefficients  to  each 
component . 

(5)  Sum  the  tangential  components  of  the  direct  and  reflected 
electric  field  at  the  observation  point. 
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D.  RESULTS  USING  THE  REFLECTION  COEFFICIENT  METHOD 


The  most  desirable  test  of  a numerical  method  is  comparison  with 
experimental  data,  and  lacking  that,  with  results  obtained  using  other 
independent  analytical  techniques.  Such  comparisons  are  presented  below 
to  validate  the  reflection  coefficient  method. 

In  figure  A-3  is  shown  the  input  resistance  of  a vertical  half-wave 
dipole  as  a function  of  height  over  a lossy  ground.  The  comparison  with 
Sommerfeld  integral  results  points  out  that  significant  discrepancies 
occur  only  when  the  dipole  is  very  close  to  the  ground  and  that  even 
when  the  tip  is  touching,  an  error  of  less  than  10  percent  is  encoun- 
tered. The  input  resistance  of  a 0. 1 X horizontal  dipole  above  a dielec- 
tric half-space  is  shown  as  a function  of  height  above  ground  in  figure 
A-4.  Comparison  with  the  scant  independent  data  indicates  close  agreement. 
A comparison  with  experimental  data  for  a purely  dielectric  ground  is 
included  in  figure  A-5.  The  plot  pertains  to  the  input  resistance  of  a 
horizontal  half-wave  dipole  as  a function  of  its  height  above  ground. 

Since  the  antenna  radiation  pattern  is  also  of  great  interest,  we  include 
in  figure  A-6  the  radiation  patterns  of  a monopole  over  a dielectric 
ground.  The  results  are  compared  to  those  presented  in  Collin  and 
Zucker  (1969)  wherein  a current  distribution  on  the  monopole  was  assumed. 
For  comparison,  a Sommerfeld  integral  calculation  is  also  provided. 

The  characteristics  of  antenna  arrays  over  lossy  ground  can  be 
easily  and  accurately  established  using  the  approximate  approach. 

Figure  A-7  contains  plots  of  the  input  resistance  and  radiation  patterns 
of  a parasitic  array  of  half-wave  dipoles.  The  discrepancies  between 
two  independent  approaches  are  noticeable  only  for  h/X  < 0.35.  The  same 
array  is  considered  with  variable  spacing  in  figure  A-8  where  a comparison 
with  the  Sommerfeld  integral  and  perfect  ground  results  is  included.  To 
emphasize  the  usefulness  of  the  method,  the  radiation  patterns  of  a log 
periodic  antenna  over  various  grounds  are  shown  in  figure  A-9-  With  a 
knowledge  of  the  free-space  or  perfect  ground  pattern,  one  must  appreciate 
that  an  accurate  prediction  of  the  pattern  over  an  arbitrary  ground 
would  indeed  be  difficult.  Furthermore,  for  the  case  shown,  the  evalua- 
tion of  the  Sommerfeld  integrals  would  be  not  only  difficult,  but  time 
consuming  and  expensive. 


A- 19 


INPUT  RESISTANCE  (ohms) 


0 0.05  0.1  0.15  0.2 
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Figure  A-**.  Input  Resistance  of  Horizontal  Dipole,  0.1X  Long,  Located 
a Distance  H Above  a Purely  Dielectric  Ground  (Dielectric 
Constant  e) 


INPUT  RESISTANCE  (ohms) 


Figure  A-6.  Radiation  Patterns  (Linear  Voltage  Scale)  of  a 
Monopole  Over  a Dielectric  Ground 
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Dipole  Length  L = 0.5X 
Separation  d =0.1X 
£x  = 10.0;  <7=  3x10" 
Frequency  = 3.0  MHz 


h/x 


a)  Input  Impedance 

■ Reflection  Coefficient 
Sommerfeld  Integral 


2 meteis  PERFECT  GROUND 

165 MHz  DIELECTRIC  GROUND 

.86  LOSSY  GROUND 

-54°  FREE  SPACE 


Figure  A-9.  Radiation  Patterns  of  a Log-Periodic 
Antenna  Over  Ground 
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METRIC  SYSTEM 


a 


BASE  UNITS: 

Quantity  Unit  SI  Symbol 


length 

metre 

m 

mass 

kilogram 

kg 

time 

second 

8 

electric  current 

ampere 

A 

thermodynamic  temperature 

kelvin 

K 

amount  of  substance 

mole 

mol 

luminous  intensity 

candela 

cd 

SUPPLEMENTARY  UNITS: 

plane  angle 

radian 

rad 

solid  angle 

steradian 

sr 

DERIVED  UNITS: 

Acceleration 

metre  per  second  squared 

activity  (of  a radioactive  source) 

disintegration  per  second 

angular  acceleration 

radian  per  second  squared 

angular  velocity 

radian  per  second 

area 

square  metre 

density 

kilogram  per  cubic  metre 

F 

electric  capacitance 

farad 

electrical  conductance 

siemens 

S 

electric  field  strength 

volt  per  metre 

H 

electric  inductance 

henry 

electric  potential  difference 

volt 

V 

electric  resistance 

ohm 

electromotive  force 

volt 

V 

energy 

joule 

1 

entropy 

joule  per  kelvin 

N 

force 

newton 

frequency 

hertz 

Hz 

illuminance 

lux 

lx 

luminance 

candela  per  square  metre 

Im 

luminous  flux 

lumen 

magnetic  field  strength 

ampere  per  metre 

Wb 

magnetic  flux 

weber 

magnetic  flux  density 

tesla 

T 

magnetomotive  force 

ampere 

A 

power 

watt 

W 

pressure 

pascal 

Pa 

quantity  of  electricity 

coulomb 

C 

quantity  of  heat 

joule 

1 

radiant  intensity 

watt  per  steradian 

specific  heat 

joule  per  kilngram-kelvin 

Pa 

stress 

past  al 

thermal  conductivity 

watt  per  metre-kelvin 

velocity 

metre  per  second 

viscosity,  dynamic 

pascal-second 

viscosity,  kinematic 

square  metre  per  second 

voltage 

volt 

V 

volume 

cubic  metre 

wavenumber 

reciprocal  metre 

work 

joule 

i 

SI  PREFIXES: 


Multiplication  Factors 


I’refix 


1 000  00(1  000  000 

in'1 

I era 

1 000  ooo  000 

10’ 

Riga 

i ooo  ooo  - 

111* 

mega 

1 ooo  = 

10’ 

kilo 

100  - 

10' 

hecto* 

10  = 

10' 

deka* 

01  = 

10-' 

deci* 

0 01  e 

10- > 

cnntr 

0 001  =-- 

10- ’ 

mill! 

0 000  001 

10-  * 

micro 

0.000  000  001 

in  * 

nano 

0.000  000  000  001  - 

10  11 

pico 

0 000  000  000  (100  001 

10-” 

fnmto 

0.000  0(H)  000  000  000  001 

10  - 

alto 

‘ To  be  avoided  where  possible 


Formula 


m/s 

(disintegration  )/s 

rad/s 

rad/s 

m 

kg/m 

A-s/V 

A/V 

V/m 

Vs/A 

W/A 

V/A 

W/A 

N-m 

|/K 

kg-m/s 

(cyde)/s 

im/m 

cd/m 

cd-sr 

A/m 

V-s 

Wb/m 

I/s 

N/m 

A-s 

N-m 

W/sr 

J/kg-K 

N/m 

W/m-K 

m/s 

Pa-s 

m/s 

W/A 

m 

( wave)/m 
N-m 


SI  Symbol 

T 

('. 

M 

k 

h 

da 

d 


m 

M 

n 


id 


H 


MISSION 

of 

Rome  Air  Development  Center 


RADC  plans  and  conducts  research , exploratory  and  advanced 
development  programs  in  command,  control,  and  communications 
( C 3)  activities , emd  in  the  C3  areas  of  information  sciences 
and  intelligence.  The  principal  technical  mission  areas 
are  communications,  electromagnetic  guidance  and  control, 
surveillance  of  ground  and  aerospace  objects,  intelligence 
data  collection  and  handling,  information  system  technology , 
ionospheric  propagation,  solid  state  sciences,  microwave 
physics  and  electronic  reliability,  maintainability  and 
compatibility . 
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